INTRODUCTION
Juniperus sabina L., commonly known as savin juniper or savin, is a morphologically variable, predominantly dioecious shrub belonging to the monophyletic section Sabina within the genus Juniperus L. (Cupressaceae) (Mao et al. 2010). It is widely distributed across mountainous regions of southern and central Europe, extending eastward into western and central Asia (Adams 2014). The species typically occupies dry, rocky habitats at mid to high elevations, where it exhibits remarkable ecological plasticity and tolerance to drought, cold, and poor soils. Despite its broad range, Juniperus sabina often forms small, scattered populations, suggesting a refugial distribution shaped by complex biogeographic and climatic histories. Previous studies have revealed weak geographical structuring in its morphology (Mazur et al. 2021), despite significant genetic differentiation across its range (Jadwiszczak et al. 2023, 2024), pointing to a long history of gene flow and relatively recent population isolation. Moreover, genome size variation and polyploidy have been documented within the species (Adams et al. 2016, 2017, 2018; Farhat et al. 2019) further complicating its taxonomic and evolutionary interpretation.
Recent taxonomic treatments recognize five infraspecific taxa within Juniperus sabina, reflecting its wide ecological amplitude and complex evolutionary history. According to the Plants of the World Online database maintained by Kew Science (2025), the accepted varieties include: var. sabina, var. balkanensis R.P.Adams et Tashev, var. arenaria (E.H.Wilson) Farjon, var. davurica (Pall.) Farjon, and var. mongolensis R.P.Adams. Of these, three varieties – davurica, mongolensis, and arenaria – are distributed across the eastern part of the species’ range in Central Asia, while two – sabina and balkanensis – are confined to Europe and western Anatolia. The European varieties are of particular interest due to their overlapping morphological traits, contrasting ploidy levels, and distinct biogeographic histories.
The most recent and taxonomically significant development concerns the recognition of Juniperus sabina var. balkanensis, a tetraploid lineage (specifically an allotetraploid) hypothesized to have originated through ancient hybridization between maternal diploid Juniperus sabina var. sabina and a paternal ancestor related to Juniperus thurifera L. (Adams et al. 2016). This variety has since been documented across the Balkan Peninsula, western Turkey, and parts of the Apennines, including newly confirmed populations in Croatia (Adams et al. 2018). Molecular studies using chloroplast DNA, nuclear microsatellites, and SNP markers have consistently supported the genetic distinctiveness of var. balkanensis (Jadwiszczak et al. 2023, 2024), despite its morphological similarity to var. sabina (Mazur et al. 2021). The near-parapatric distribution of the two varieties, coupled with their contrasting ecological preferences – particularly in relation to soil potassium content and drought tolerance – suggest a complex evolutionary trajectory shaped by hybridization, polyploidization, and climatic filtering. Var. sabina occupies cooler and wetter environments, whereas var. balkanensis is adapted to warmer, drier habitats with potassium-rich soils (Jadwiszczak et al. 2023, 2024).
In Croatia, Juniperus sabina var. balkanensis forms dis-tinctive shrubland communities classified within the EUNIS habitat system as F2.2 – evergreen alpine and subalpine heath and scrub. These habitats are typically found above 900 m a.s.l. in the Dinaric Alps (Figure 1), occupying xerothermic karst terrains with shallow, skeletal soils developed on limestone and dolomite. The Croatian populations of Juniperus sabina var. balkanen-sis are characterized by dense, prostrate growth forms that inhibit the development of continuous grass cover, thereby creating patchy microhabitats within mosaics of dry grasslands and rocky outcrops. The vegetation includes a suite of xerophilous and thermophilus species such as Juniperus communis var. saxatilis Pall., Amelanch-ier ovalis (Willd.) Borkh., Rhamnus saxatilis Jacq., and Arctostaphylos uva-ursi (L.) Spreng., alongside grassland elements like Sesleria juncifolia Suffren, Festuca spp., and Stipa spp. These shrublands represent a transitional suc-cessional stage toward Fagus sylvatica L. or Pinus nigra J.F. Arnold forests, with Juniperus sabina var. balkan-ensis playing an important role in soil stabilization and organic matter accumulation through its extensive root system and needle litter (Topić and Vukelić 2009).
Although the presence of the tetraploid Juniperus sabina var. balkanensis in Croatia has been confirmed (Adams et al. 2018), its morphological and anatomical variability along the Dinaric range remains poorly understood. Previous studies have focused exclusively on a single population from the Velebit Mountain (Mazur et al. 2021), leaving a significant gap in our knowledge of the species’ intraspecific diversity in this region. In the present study, we examined three natural populations – two from different localities on the Velebit Mountain and one from the Biokovo Mountain – representing the southernmost and ecologically contrasting segments of the Croatian Dinarides. In contrast to the widespread Juniperus communis var. saxatilis, which dominates subalpine scrublands across the Dinaric karst, Juniperus sabina var. balkanensis is relatively rare and spatially restricted. Its occurrence in Croatia is limited to isolated patches within high-altitude, xerothermic habitats, making it a valuable element of regional biodiversity. The Dinaric Alps, with their pronounced ecological gradients and floristic heterogeneity, thus provide a unique setting for investigating population-level variation in relict and alpine taxa.
This study aimed to assess the morphological and anatomical variability of vegetative traits (scale leaves and shoots), as well as the morphological diversity of cones within and among populations. Special attention was given to differences between male and female shrubs based on vegetative traits. Furthermore, we explored the extent to which population structure is influenced by geographic position, elevation, and other environmental factors, with the goal of identifying patterns of local adaptation and potential drivers of phenotypic divergence.
Figure 1 Habitat of the Biokovo population, showing scattered shrubs of Juniperus sabina var. balkanensis among dominant shrubs of Juniperus communis var. saxatilis.
MATERIALS AND METHODS
Study area and plant material
In the spring of 2025, plant material of Juniperus sabina var. balkanensis was collected for anatomical and morphometric research from three natural populations in Croatia (Figure 2, Table 1): Veliki Alan and Mali Alan on the Velebit Mountain, and one on the Biokovo Mountain. A total of 35 individual shrubs were sampled – 13 from Veliki Alan, 12 from Mali Alan, and 10 from Biokovo – comprising 16 female and 19 male plants (Table 1). Each shrub was georeferenced using GPS to ensure precise spatial documentation and support future environmental interpretation. From every shrub, a single branch approximately 20–30 cm in length was clipped to examine the anatomical structure and morphometry of scale-leaves. On female plants, berry-like seed cones were also collected, with around 10 cones taken per shrub. All samples were immediately placed in ziplock plastic bags or small containers in the field to preserve freshness and prevent desiccation. Each sample was labelled with a unique identifier to ensure accurate tracking. The collected material was then transported to the Laboratory for Morphological Research at the Department of Forest Genetics, Dendrology and Botany, University of Zagreb Faculty of Forestry and Wood Technology, for further analysis. Prior to analysis, all plant material was stored under refrigeration to maintain sample integrity. Additionally, representative specimens from each locality were prepared as herbarium vouchers and deposited in the DEND Herbarium at the University of Zagreb Faculty of Forestry and Wood Technology.
Figure 2 Geographic distribution of three sampled populations of Juniperus sabina var. balkanensis in the Croatian Dinaric
Alps, with the overall range of the taxon in Croatia shown in blue.
Table 1 Populations of Juniperus sabina var. balkanensis from the Croatian Dinaric Alps, with geographic coordinates, elevation, number of analysed shrubs (total, male and female), and bioclimatic variables (BIO1, BIO3, BIO4, BIO12, BIO15, SOLAR) used in the Redundancy Analysis (RDA). Bioclimatic variables were calculated as population means based on individual sampling coordinates, while elevation is presented as mean values with ranges. Although BIO12 (annual precipitation) is shown in the table to facilitate interpretation, it was excluded from Redundancy Rnalysis (RDA) due to collinearity. Bioclimatic variables: BIO1 (mean annual temperature), BIO3 (isothermality), BIO4 (temperature seasonality), BIO12 (annual precipitation), BIO15 (precipitation seasonality), and SOLAR (annual solar radiation).
Morphological and anatomical characterization of shoots, scale leaves, and cones
Morphological analysis of scale leaves and shoots was performed using both stereo and compound microscopy. Macrophotographs of freshly collected material were obtained with a stereo microscope (Olympus SZX12) equipped with a digital camera (Olympus EP50), while microphotographs of leaf cross-sections were taken with a compound microscope (Olympus BX41) fitted with a digital camera (Olympus LC35). Vegetative morphometric traits – including the number of lateral (NLS) and facial (NFS) scale leaf pairs per 1 cm of shoot length, facial leaf length (FSL), and overall shoot width (SW) – were measured at a standardized distance of 0.5 cm from the shoot apex. These measurements were performed using the integrated On-Screen Display (OSD) software of the Olympus EP50 camera. For anatomical analysis, epidermal thickness (ET) was quantified using PRECiV Capture software (version 1.2, Build 27732) (PRECiVTM 2021), based on calibrated microphotographs of transverse leaf sections.
Cone dimensions were measured using a digital calliper (Alpha Tools®, Bahag AG, Germany, ±0.1 mm). For each cone, three measurements were taken: cone length (CL), cone width (CW), and cone thickness (CT). After these morphometric measurements, the total number of seeds contained in each cone (NS) was determined.
Statistical analysis
Morphological variability of shoots, scale leaves and seed cones was assessed using descriptive statistical parameters. For each trait, the arithmetic mean, standard deviation, minimum and maximum values, and the coefficient of variation (CV, %) were calculated. For scale leaf traits these metrics were computed separately for male and female individuals, as well as for each population and for the total dataset.
Morphometric traits of vegetative shoots and cones were analysed using linear mixed-effects models in R (R Core Team 2025). The analysis was performed using the lme4 package (Bates et al. 2025), which allows specification of both fixed and random effects. For shoot-related traits, the model included fixed effects of sex and population, as well as their interaction, while shrub nested within population was treated as a random effect. For cone-related traits, fixed effects included population, and shrub nested within population was modelled as a random effect. Model fitting was conducted using the lmer() function, and the significance of fixed effects was assessed via Type III ANOVA using the car package (Fox and Weisberg 2019).
Principal Component Analysis (PCA) was conducted to examine patterns of variation among shrubs and populations based on morphometric traits. PCA was performed separately for vegetative and cone morphological traits using the FactoMineR (Lê et al. 2008) and Factoextra (Kassambara and Mundt 2020) packages in R Statistical Software (v4.4.0; R Core Team 2025). Prior to analysis, all trait variables were standardized (z-transformed) to ensure comparability and reduce scale effects. For vegetative traits, PCA was based on five variables, while cone traits were analysed using four. Eigenvalues, explained variance, and variable loadings were extracted to assess component significance and trait contributions. Components with eigenvalues greater than one were retained for interpretation. Population-level patterns were assessed based on individual scores along the first two principal components, and visualized using biplots that simultaneously display trait orientation and individual clustering.
To assess the influence of environmental and geographic factors on morphological variation in Juniperus sabina var. balkanensis, we performed Redundancy Analysis (RDA) using the vegan package (v2.6-4; Oksanen et al. 2022) within R Statistical Software (v4.4.0; R Core Team 2025). Morphological traits were divided into two sets: vegetative traits and cone traits, each analysed separately using constrained ordination. Prior to analysis, all morphological variables were Hellinger-transformed to reduce scale dependency and meet assumptions of linear ordination. Traits with high intercorrelation (r > 0.85) were excluded to avoid redundancy, retaining only those with clear population-level differentiation. Environmental predictors included bioclimatic variables (BIO1–BIO19), solar radiation, and altitude. Bioclimatic data were retrieved from WorldClim v2.1 (Fick and Hijmans 2017) at 30 arc-second resolution, while altitude was obtained via GPS measurements. To reduce multicollinearity, we performed a correlation analysis among environmental variables and excluded all variables with pairwise Pearson correlation coefficients exceeding ±0.9. Final predictor sets were selected separately for each trait group based on ecological relevance and statistical independence. All predictors were standardized (z-transformed) prior to ordination. Geographic predictors included latitude and longitude. To partition the variation attributable to environment and geography, we constructed three RDA models for each trait group: (1) a full RDA model including both environmental and geographic predictors; (2) a partial RDA isolating environmental effects conditioned on geography (environment | geography); (3) a partial RDA isolating geographic effects conditioned on environment (geography | environment). Model significance was assessed using permutation tests (999 iterations), and the adjusted R² values were used to evaluate explanatory power.
Ordination results from both PCA and RDA analyses were visualized using the ggplot2 package (Wickham 2016) in R Statistical Software (v4.4.0), with enhancements provided by ggrepel (Slowikowski 2024) for label clarity and ggsci (Xiao 2024) for consistent colour schemes. All plots were standardized for axis scaling and annotated to highlight population-level patterns and trait contributions. Only the first two axes were displayed in each ordination, as they captured the majority of explained variation.
RESULTS
Variation in vegetative traits
Across all individuals, scale leaf traits exhibited low to moderate variability (Table 2). The number of facial scale leaf pairs (NFS) ranged from 5 to 15, with a mean of 7.40 ± 1.15, while lateral scale leaf pairs (NLS) ranged from 4.5 to 15, averaging 7.26 ± 1.16 per centimetre of shoot length. Shoot width (SW) was relatively stable (1.32 ± 0.12 mm, range: 1.03–1.72 mm), whereas facial leaf length (FSL) showed greater variation (1.65 ± 0.28 mm, range: 1.05–2.84 mm). Epidermal thickness (ET) ranged from 11.00 to 22.18 μm, with a mean of 16.79 ± 1.62 μm. Coefficients of variation (CV) below 10% were observed for shoot width (8.73%) and epidermal thickness (9.67%), while higher CVs were recorded for facial scale leaf number (15.53%), lateral scale leaf number (16.00%), and scale leaf length (16.68%).
Morphological traits of scale leaves varied across the three populations – Veliki Alan, Mali Alan, and Biokovo – with notable differences in both size and number (Table 2). The Biokovo population exhibited the highest mean values for both facial (8.09) and lateral (7.92) scale leaf pairs, with broader ranges (facial: 6–15; lateral: 6–15) and slightly higher variability in lateral leaf count (CV = 15.30%). In contrast, Mali Alan had the lowest mean values (facial: 7.05, lateral: 6.87), though variability remained comparable (CV ≈ 14.3–14.6%). Shoot width was consistent across populations, ranging from 1.31 to 1.34 mm, with low CVs (<10%). The shortest facial scale leaves were recorded in the Biokovo population (mean = 1.50 mm, CV = 11.35%), while Mali Alan had the longest scale leaves (mean = 1.73 mm, CV = 16.11%). Veliki Alan showed similar values to Mali Alan (mean = 1.70 mm, CV = 17.24%), indicating a shared trend in leaf elongation. Epidermal thickness was highest in Biokovo (mean = 17.21 μm, range: 11.00–22.18 μm), with slightly greater variability (CV = 11.65%) compared to Veliki Alan (mean = 16.46 μm, CV = 8.52%) and Mali Alan (mean = 16.78 μm, CV = 8.30%).
Variation in seed cone traits
Among the four cone traits analysed across studied populations (Table 3), the number of seeds per cone (NS) emerged as the most variable trait (CV = 42.65%). In contrast, cone length (CL) and cone width (CW) showed relatively low variability (CV = 10.93% and 11.53%, respectively), suggesting morphological stability in these dimensions. Cone thickness (CT) displayed intermediate variability (CV = 16.22%).
Among the four cone traits analysed across Juniperus sabina var. balkanensis populations (Table 3), cone length (CL) showed the most pronounced differences in mean values. Biokovo had significantly shorter cones (mean = 5.47 mm) compared to Veliki Alan (6.16 mm) and Mali Alan (6.41 mm). In contrast, cone width (CW) remained relatively consistent across populations, ranging from 6.11 to 6.36 mm, and cone thickness (CT) from 5.50 to 5.82 mm. The number of seeds per cone (NS) was highest in Biokovo (mean = 2.42), followed by Mali Alan (2.12) and Veliki Alan (1.84), suggesting slight reproductive differentiation among sites.
Sexual dimorphism
Sex-based differences were observed across several measured traits (Table 2). The number of facial scale leaf pairs (NFS) was consistently higher in female shrubs across all populations. For example, in the Biokovo population, female shrubs averaged 8.52 pairs compared to 7.66 in male shrubs, while in the Veliki Alan population the difference was 7.59 vs. 6.76. A similar pattern was found for lateral scale leaves (NLS), with mean values in female shrubs ranging from 7.11 to 8.42, and in male shrubs from 6.69 to 7.42.
Shoot width (SW) showed minimal variation between sexes, with nearly identical mean values (1.31–1.34 mm) and low coefficients of variation (CV < 10%) in both groups. In contrast, facial scale leaf length (FSL) was consistently greater in male shrubs. In the Mali Alan population, mean FSL was 1.77 mm in male shrubs compared to 1.69 mm in female shrubs, whereas in Veliki Alan the difference was even more pronounced (1.82 mm vs. 1.58 mm).
Table 2 Descriptive statistics of vegetative traits in Juniperus sabina var. balkanensis across populations and sexes. Descriptive parameters: –x – arithmetic mean; Min – minimum value observed within the sample; Max – maximum value observed within the sample; SD – standard deviation; CV (%) – coefficient of variation. Acronyms of vegetative morphometric traits: NFS – number of facial scale pairs; NLS – number of lateral scale pairs; SW – shoot width; FSL – facial scale length; ET – epidermis thickness.
Table 3 Descriptive statistics for analysed cone traits for Juniperus sabina var. balkanensis. Descriptive parameters: –x – arithmetic mean; Min – minimum value observed within the sample; Max – maximum value observed within the sample; SD – standard deviation; CV (%) – coefficient of variation. Acronyms of cone morphometric traits: CL – cone length; CW – cone width; CT – cone thickness; NS – number of seeds per cone.
Epidermal thickness (ET) was slightly higher in female shrubs across all populations. In the Biokovo population, female shrubs reached a mean of 17.27 μm, compared to 17.15 μm in male shrubs, with similar trends observed in Veliki Alan (16.84 μm vs. 16.09 μm). In the Mali Alan population, however, epidermal thickness was nearly identical between female and male shrubs (16.78 μm vs. 16.77 μm).
Results of the General Linear Model (GLM) for veg-etative and cone traits
Populations and sexes differed significantly in three traits (Table 4): the number of facial scale leaf pairs (NFS), the number of lateral scale leaf pairs (NLS), and the facial scale leaf length (FSL). Population effects explained 16.94%, 15.12%, and 14.45% of the total variance for these traits, respectively, while sex contributed 14.52%, 13.37%, and 6.59%. In contrast, shoot width (SW) and epidermal thickness (ET) showed no significant differences between populations or sexes, with population-level effects explaining only 0–3.15% and sex effects below 0.40%. These traits appear morphologically stable across broader spatial and habitat gradients.
The interaction between population and sex was non-significant in all cases (Table 4), contributing 0% to the total variance, suggesting that sex-related differences were consistent across populations and not context-dependent.
Shrubs nested within populations accounted for a substantial portion of the variance across all traits (Table 4), ranging from 15.02% (ET – epidermal thickness) to 28.42% (NLS – the number of lateral scale leaf pairs), thus highlighting pronounced individual-level variability within populations. Residual error explained the remaining variation, ranging from 43.09% to 81.44%, depending on the trait, likely reflecting both measurement noise and unaccounted biological complexity (e.g., within-shrub variability).
Table 4 Results of the General Linear Model (GLM) for vegetative traits. Acronyms of vegetative morphometric traits: NFS – number of facial scale pairs; NLS – number of lateral scale pairs; SW – shoot width; FSL – facial scale length; ET – epidermis thickness. df – degrees of freedom; F – F-test in GLM; p – significance of GLM’s F-test. *** significant at p < 0.001; ** significant at 0.001 < p < 0.01; * significant at 0.01 < p < 0.05; ns depicts non-significant values (p > 0.05).
Significant variation in cone length (CL) was observed among the studied Juniperus sabina var. balkanensis populations (Table 5), accounting for 43.68% of the total variance, while within-population variation was considerably lower (8.50%). In contrast, cone width (CW) and cone thickness (CT) showed no significant inter-population differences, with most of the variation attributed to within-population variability and residual error. The number of seeds per cone (NS) exhibited relatively low within-population variation (10.91%), whereas differences among populations were not statistically significant. These findings indicate that cone length (CL) is the most spatially and environmentally structured trait, while cone width (CW) and thickness (CT) are primarily influenced by individual-level variation within populations.
Principal Component Analysis (PCA)
Principal Component Analysis (PCA) based on five vegetative traits revealed that the first two principal components had eigenvalues greater than one and together explained 80.3% of the total variance. The first principal component (PC1) accounted for 58.0% of the variation and was strongly positively correlated with the number of facial scale leaves pairs (NFS, r = 0.974) and lateral scale leaves pairs (NLS, r = 0.944), while showing a strong negative correlation with facial scale leaf length (FSL, r = -0.907). The second principal component (PC2), which explained an additional 22.3% of the variance, was primarily associated with shoot width (SW, r = -0.831) and epidermal thickness (ET, r = -0.594), both showing moderate to strong negative loadings. The remaining components (PC3–PC5) had eigenvalues below one and together explained the remaining 19.7% of the total variance (PC3 = 15.1%, PC4 = 4.3%, PC5 = 0.3%), but were not considered informative for interpretation.
Table 5 Results of the General Linear Model (GLM) for cone traits. Acronyms of cone morphometric traits: CL – cone length; CW – cone width; CT – cone thickness; NS – number of seeds per cone. df – degrees of freedom; F – F-test in GLM; p – significance of GLM’s F-test. *** significant at p < 0.001; ** significant at 0.001 < p < 0.01; * significant at 0.01 < p < 0.05; ns depicts non-significant values (p > 0.05).
Despite partial overlap among Juniperus sabina var. balkanensis shrubs, the PCA biplot revealed a clear spatial separation of populations (Figure 3A). Individuals from the Biokovo population clustered predominantly on the right side of the diagram, characterized by higher numbers of facial (NFS) and lateral (NLS) scale leaf pairs. In contrast, shrubs from Mali Alan and Veliki Alan (Velebit Mountain) populations grouped on the left side, distinguished by longer facial scale leaves (FSL).
Similar to the PCA based on vegetative traits, the analysis of cone morphology also revealed structured variation among populations (Figure 3B). The first two principal components had eigenvalues greater than one and together explained 81.7% of the total variance. PC1 accounted for 50.5% and was strongly negatively correlated with cone thickness (CT, r = –0.864), cone width (CW, r = –0.800), and the number of seeds per cone (NS, r = –0.756), indicating that individuals with higher scores on this axis tend to have smaller cones with fewer seeds. PC2 explained an additional 31.2% of the variance and was strongly negatively associated with cone length (CL, r = –0.934), while also showing a moderate positive correlation with seed number (NS, r = 0.578). The remaining components (PC3 and PC4) had eigenvalues below one and together accounted for the remaining 18.3% of the variance, indicating that most of the morphological variation in cone traits is captured by the first two axes. Unlike the PCA with vegetative traits, where population separation was more gradual, the biplot based on cone traits revealed a pronounced divergence along the second principal component. Shrubs from the Biokovo population formed a clearly distinct cluster on the upper side of the diagram, characterized by shorter cones with a larger number of seeds, whereas individuals from Mali Alan and Veliki Alan (the Velebit Mountain) grouped together on the lower side, defined by longer cones with fewer seeds per cone.
Redundancy Analysis (RDA)
To examine how environmental and spatial gradients influence vegetative morphology in Juniperus sabina var. balkanensis, we performed Redundancy Analysis (RDA) using eight explanatory variables: six ecological (BIO1 – mean annual temperature, BIO3 – temperature isothermality, BIO4 – temperature seasonality, BIO15 – precipitation seasonality, solar radiation and altitude) and two geographical (latitude and longitude). Of the five vegetative traits initially considered, four were retained for analysis due to high collinearity between the number of lateral (NLS) and facial (NFS) scale leaf pairs. The full model approached statistical significance (F = 2.009, p = 0.058) and explained 38.2% of the variance observed in vegetative traits (Table 6). The first constrained axis (RDA1) accounted for 85.9% of the explained variation, while the second axis (RDA2) contributed an additional 13.5% (Figure 3C). The remaining axes explained negligible amounts and were not further considered. Partial RDA models revealed that environmental variables (pRDAenv) accounted for a larger portion of the variation (R² = 0.206) than geographical variables (pRDAgeo, R² = 0.011), which were not significant. Several environmental and geographical predictors showed statistically significant associations with vegetative morphological variation, as indicated by permutation tests (p < 0.05). These included temperature isothermality (BIO3), temperature seasonality (BIO4), solar radiation (SOLAR), latitude, and longitude. Marginally significant effects were observed for precipitation seasonality (BIO15) and altitude (p < 0.1), while mean annual temperature (BIO1) was not significant. In the ordination space (Figure 3C), most of these significant predictors were strongly aligned with the first constrained axis (RDA1), particularly BIO3, BIO4, and SOLAR, which exhibited high negative loadings. Latitude and longitude also showed strong projections along RDA1, suggesting that both environmental and spatial gradients contribute to the structuring of vegetative traits. RDA2 captured more subtle variation, with weaker associations to the included predictors.
In the second model, we examined how cone morphol-ogy responds to ecological and spatial gradients (Figure 3D, Table 7). Te model included five explanatory varia-bles: mean annual temperature (BIO1), temperature iso-thermality (BIO3), altitude, latitude, and longitude. The full RDA model was statistically significant (F = 4.1784, p = 0.001), explaining 67.6% of the variability observed in cone traits (adjusted R² = 0.514). Variance partition-ing using partial RDA revealed that environmental vari-ation independent of geography accounted for 24.3% of the trait variance (adjusted R² = 0.169, p = 0.001), while the unique effect of geography – after controlling for en-vironmental factors – explained only 13.1% and was not statistically significant (adjusted R² = 0.083, p = 0.14). Among the predictors, the environmental variable tem-perature isothermality (BIO3), along with latitude and longitude, exhibited strong and statistically significant associations with cone traits (p < 0.01). These predic-tors were predominantly aligned with the first RDA axis, which captured 92.2% of the explained variation. In con-trast, RDA2 accounted for only 6.4% and did not show strong associations with any individual predictor.
Figure 3 Multivariate ordination analyses (PCA and RDA) of vegetative and cone traits in Juniperus sabina var. balkanensis populations from the Croatian Di-naric Alps. (A–B) Principal Component Analyses (PCA) showing differentiation among populations based on vegetative (A) and cone traits (B). Individuals are represented as coloured points grouped by population, with arrows indicating the direction and magnitude of trait loadings. (C–D) Redundancy Analyses (RDA) illustrating population-level variation in vegetative (C) and cone traits (D) in relation to environmental and geographic gradients. Circles denote individuals grouped by population, squares represent morphological traits, and vectors show the strength and orientation of environmental and geographical predictors. Acronyms for environmental variables are defined in Table 1; acronyms for morphological and anatomical traits are listed in Tables 2 and 3.
Table 6 Performance metrics for full and partial RDA models assessing the effects of environmental and geographic predictors on vegetative morphological variation in Juniperus sabina var. balkanensis. Adjusted R², F-values, and p-values are based on permutation tests with 999 iterations.
Table 7 Performance metrics for full and partial RDA models assessing the effects of environmental and geographic predictors on cone morphological variation in Juniperus sabina var. balkanensis. Adjusted R², F-values, and p-values are based on permutation tests with 999 iterations.
DISCUSSION
Although our data fall within the broader morphological variation reported for Juniperus sabina in dendrological literature (Krüssmann 1972, Vidaković 1991, Farjon 2005, Idžojtić 2013, Adams 2014), direct comparisons with earlier studies remain challenging due to methodological inconsistencies – particularly in how scale leaves and shoot segments are defined and measured. As a result, traits related to vegetative morphology are difficult to compare reliably across studies, whereas cone dimensions – being less sensitive to morpho-anatomical interpretation – allow for more meaningful comparisons.
For instance, Amaral Franco (1964) and Mazur et al. (2021) reported notably smaller shoot width values (0.6–0.8 mm and 0.9 mm, respectively), both of which fall below the average observed in our populations (1.3 mm). These discrepancies likely stem from differences in sampling protocols, measurement criteria, and the type of material analysed. Mazur et al. (2021) explicitly stated that their material was dried without pressing prior to measurement, whereas our study relied on freshly collected specimens. Drying can cause shrinkage and deformation of vegetative structures, especially scale leaves and thin shoots, which likely contributed to the smaller values reported in their study. In addition, discrepancies may also arise from differences in the shoot segment selected for measurement. A similar divergence was observed in leaf counts: Mazur et al. (2021) reported approximately 18 scale leaves per 0.5 cm of shoot, while our counts reached 30 scale leaves per centimetre – equating to roughly 15 per 0.5 cm. Finally, the average length of facial scale leaves in our study was 1.6 mm, which fits well within the 1–3 mm range reported by Krüssmann (1972). He emphasized that leaf scale size varies depending on shoot type – being shorter on lateral branchlets and longer on apical axes – which further supports the consistency of our sampling, based on facial scale leaves measured on short lateral branchlets. Altogether, these comparisons highlight the importance of clearly defined sampling criteria in morphometric studies. Given the developmental plasticity and structural variability of vegetative traits in juniper species with scale leaves (Mazur et al. 2003, Douaihy et al. 2012, Boratyński et al. 2013), standardized definitions, consistent measurement protocols, and careful consideration of material type (fresh vs. dried) are essential for ensuring reproducibility and meaningful comparisons across studies.
Cone traits showed general agreement with published data (Krüssmann 1972, Vidaković 1991, Farjon 2005, Idžojtić 2013, Adams 2014), particularly with the values reported for Juniperus sabina in Mazur et al. (2021). In our populations, mean cone length (CL) was 6.0 mm, cone width (CW) averaged 6.1 mm, and cone thickness (CT) reached 5.8 mm. These values are slightly lower than those reported for var. balkanensis (CL = 6.5 mm, CW = 6.8 mm, CT = 6.0 mm), but closely resemble those recorded for var. sabina (CL = 6.1 mm, CW = 6.0 mm, CT = 5.3 mm), indicating partial overlap in cone morphology between the two varieties. The average number of seeds per cone in our samples was 2.1, nearly identical to the mean of 2.2 reported by Mazur et al. (2021), and within the species-level range described by Krüssmann (1–3 seeds). In our case, seed number ranged from 1 to 4, suggesting slightly broader variability. Likewise, Krüssmann (1972) noted that cone length typically falls between 5 and 7 mm, which aligns well with our observed values.
While our measurements generally fall within the ranges reported in literature (Krüssmann 1972, Vidaković 1991, Farjon 2005, Idžojtić 2013, Adams 2014, Mazur et al. 2021), the morphological distinction between Juniperus sabina var. balkanensis and var. sabina remains weakly defined. As emphasized by Mazur et al. (2021), most traits showed statistically significant differences between varieties, yet the overall morphological structure was only weakly expressed. This supports the interpretation of var. balkanensis as morphologically cryptic, with its recognition based primarily on genetic evidence (Jadwiszczak et al. 2023, 2024), particularly plastid DNA indicative of ancient hybridization events, as first described by Adams et al. (2016). Our findings reinforce this view: morphological traits in var. balkanensis populations do not provide clear diagnostic separation from var. sabina described in the literature (Adams et al. 2016). Taken together, these results highlight the limitations of morphology alone in resolving intraspecific variation of Juniperus sabina and underscore the need for integrative approaches that combine morphometric, ecological, and molecular data.
Sex-based differences were observed in several morphometric traits, most notably in the number and length of scale leaves. Female shrubs consistently produced more facial and lateral scale leaves per centimetre of shoot, whereas male shrubs exhibited slightly longer facial scale leaves. Epidermal thickness was also marginally greater in female shrubs across populations. Shoot width showed no variation between sexes, and none of the traits exhibited a significant population × sex interaction, indicating that sexual dimorphism was consistent across populations and not context-dependent. Nonetheless, with a sample size of 35 shrubs from only three populations, these findings should be interpreted with caution regarding broader species-level patterns. Although female plants produced more scale leaves per centimetre of shoot, this does not necessarily imply a greater vegetative investment overall, as total shoot volume and biomass were not assessed. In gymnosperms, female plants often develop larger or more numerous leaves to support reproductive demands (Gauquelin et al. 2002, Nowak-Dyjeta et al. 2017, Stefanović et al. 2017, Vidaković et al. 2024), particularly cone maturation, which requires substantial assimilative resources. Our results partially align with this pattern, as female shrubs exhibited a higher number of scale leaves and slightly thicker epidermis. However, interpretation is further complicated by the ambiguous sexual expression of Juniperus sabina, a species known to exhibit both monoecious and dioecious individuals, and by the fact that male plants in this study tended to produce slightly longer facial scale leaves. Field observations revealed individuals bearing exclusively male cones, exclusively female cones, and some bearing both. All plants from which the seed cones were collected for morphometric analysis were classified as female, although some exhibited both male and female cones. This plasticity in sex expression may have influenced the final results. In this context, the reduced leaf length observed in female plants could reflect the energetic cost of producing both male and female reproductive structures. Taken together, our results suggest that sexual dimorphism in Juniperus sabina var. balkanensis is present but nuanced, shaped by reproductive demands, ecological interactions, and potentially flexible sex expression. Future studies with broader sampling and integrated physiological measurements are needed to clarify the functional significance of these patterns.
Across all populations, both our study and that by Mazur et al. (2021) revealed low to moderate coefficients of variation, indicating relatively stable trait expression and limited within-population morphological variability. Nonetheless, statistical analyses confirmed that intra-population variability was significant for all examined traits. These differences may reflect the underlying genetic heterogeneity among shrubs (Jadwiszczak et al. 2023, 2024), as well as individual-level plasticity shaped by microhabitat conditions such as soil composition, light exposure, and moisture availability (Douaihy et al. 2012, Boratyński et al. 2013, Vidaković et al. 2024). Although our study was based on morphometric data, the observed patterns of trait stability and low to moderate intra-population differentiation are broadly consistent with the genetic diversity revealed by SilicoDArT and SNP analyses (Jadwiszczak et al. 2023), as well as microsatellite data from populations from Europe and Asia (Jadwiszczak et al. 2024). This moderate variability at the intra-population level may be shaped by processes such as inbreeding, clonal propagation, and limited gene flow between fragmented populations. These factors can constrain genetic exchange and reduce overall diversity, particularly in small or geographically isolated stands.
As in the study by Mazur et al. (2021), our results confirmed that populations of Juniperus sabina var. balkanensis can be differentiated based on both vegetative and cone traits. In our dataset, statistically significant differences were observed for three out of five vegetative traits and for one out of four cone traits. Among all analysed traits, cone length proved to be the most informative for population-level differentiation. Regarding vegetative traits, the number of lateral and facial scale leaf pairs, as well as the length of facial scale leaves, showed consistent variation among populations. In contrast, shoot width and epiderma thickness did not differ significantly, suggesting that these traits may be less sensitive to local ecological variation or more developmentally constrained. Interestingly, these findings partially diverge from those reported by Mazur et al. (2021), who identified cone thickness, shoot width, and seed length as the most discriminative traits. The discrepancy between studies likely reflects differences in the sampling area. While Mazur et al. (2021) included a broader range of populations across southeastern Europe, our study was limited to three populations within the Dinaric region of Croatia.
The clear morphological differentiation of the Biokovo population from those on the Velebit Mountain, evident in both PCA and RDA analyses, highlights ecologically and geographically structured variation within Juniperus sabina var. balkanensis in the studied area. Remarkably, such divergence occurs across a relatively small geographic scale: the Biokovo Mountain lies only about 150 km in a straight line from the Velebit Mountain, yet the two ranges are separated by extensive stretches of unsuitable habitat for the studied species. This barrier likely restricts gene flow and reinforces isolation, providing a framework for understanding the observed differentiation. Our results showed that in both cases – whether vegetative organs or cones were considered – the full RDA models were significant, confirming that trait variation is best explained when environmental and spatial gradients are combined. Geographically isolated and exposed to higher solar radiation and reduced precipitation compared to northern ranges (Table 1), the Biokovo population exhibited a greater number of scale leaves and shorter cones with more seeds, whereas the Velebit populations had longer facial scale leaves in lower numbers and longer cones with fewer seeds.
These contrasting trait patterns, shaped by ecological gradients and reinforced by geographic isolation, may promote long-term divergence through local adaptation and phenotypic plasticity. Previous genetic studies have already demonstrated that Juniperus sabina populations are significantly differentiated and that gene flow across mountain ranges is strongly limited (Jadwiszczak et al. 2023, 2024), resulting in an isolation-by-distance pattern. Our findings of morphological divergence and geo-ecological structuring are therefore consistent with this broader genetic evidence, reinforcing the interpretation that differentiation in Juniperus sabina var. balkanensis results from the restricted gene flow imposed by mountain ranges.
CONCLUSIONS
This study provides a detailed assessment of vegetative and cone trait variation in Juniperus sabina var. balkanensis across three natural populations in the Croatian Dinaric Alps. Overall, the studied traits exhibited low to moderate variability within populations, consistent with previous genetic studies suggesting limited gene flow and clonal propagation.
Facial and lateral scale leaves differed among populations in both number and size. The Biokovo population showed the highest mean values and broader ranges for leaf scale number, whereas populations from the Velebit range (Veliki Alan and Mali Alan) were characterized by longer facial scale leaves. Shoot width and epidermal thickness remained relatively stable across sites. Within-population variability was most pronounced in leaf number and length, with individual shrubs contributing substantially to overall morphological diversity.
Cone traits exhibited the clearest and most significant differentiation in cone length, with Biokovo shrubs producing shorter cones on average, while populations from the Velebit range developed longer ones. By contrast, cone width and thickness varied predominantly at the individual level and did not show consistent spatial structuring. Seed number per cone contributed further to the observed variability: Biokovo cones, although shorter, contained more seeds, whereas Velebit cones were longer but held fewer seeds.
Sexual dimorphism was evident in several vegetative traits: female plants consistently produced more numerous and slightly thicker scale leaves, while males developed longer facial leaves. These patterns were uniform across populations, suggesting intrinsic sex-linked developmental differences rather than site-specific effects.
Multivariate analyses provided further insight into population variability. PCA highlighted the traits most responsible for population differentiation – leaf number and length among vegetative traits, and cone size and seed number among generative traits – while RDA confirmed that these differences are not random but structured by environmental and spatial gradients.
Taken together, these findings underscore the impor-tance of variation at multiple levels – among individu-al shrubs within populations, across populations, and between sexes – alongside the combined influence of ecological and spatial gradients in shaping trait expres-sion within this morphologically cryptic lineage. The observed patterns highlight the need for integrative ap-proaches combining morphometric and ecological data to better understand the evolutionary dynamics of Juni-perus sabina var. balkanensis and to inform conservation strategies for fragmented mountain shrublands.