Weak Geographical Structure of Juniperus sabina (Cupressaceae) Morphology despite Its Discontinuous Range and Genetic Differentiation

In Europe, Juniperus sabina L. is a mountainous, rare species that creates small, scattered populations, suggesting their refugial nature. Recently, a new variety of this juniper, J. sabina var. balkanensis R. P. Adams et A. N. Tashev was described based on genetic studies. We expected morphological differentiation among isolated parts of the species range and between varieties, as was the case with other Mediterranean junipers. Cones, seeds and fragments of shoots from a total of 506 individuals were collected from 24 populations in Europe and for comparisons from three populations from Tian Shan. Almost all of the 16 analysed features signiﬁcantly differentiated among populations and geographical regions as well as between the varieties, although most groups differed from others only in terms of a single feature. The thickness of cones, the width of shoots and the length of seeds were the most important features for differentiation. The geographical structure of the variation of J. sabina was weak, and comparative populations from Tian Shan were clustered with European populations, similar to the ﬁndings of a previous study on essential oils. We found slightly different patterns of variation of the two varieties of the species. The little intra-species differentiation could be the result of the long period of contact between nowadays distinct populations and their relatively late separation in the early Holocene.


Introduction
According to the latest classification, Savin juniper (Juniperus sabina L. s.l.) is divided into two species, namely the western J. sabina L. (=J. sabina var. sabina sensu Farjon [1]) and the eastern J. davurica Pall., which is further differentiated into three varieties [2]. The species belongs to the monophyletic Sabina section, the largest one within the genus Juniperus L. [3]. The range of typical J. sabina specimens stretches along 6000 km, from the westernmost localities on the Iberian Peninsula to the easternmost ones in Mongolia [4]. In the western part of its distribution range, J. sabina grows in several isolated mountainous areas of the Mediterranean region and Central Europe. In Spain, the species can be found in the Iberian System, the Cantabrian Mountains, the Betic System and, albeit more scarcely, in the Pyrenees [5,6], in the Balkans-in the Dynarides, Šar Planina, the Rodopes, Rila Mountains and Kožuf [5,[7][8][9]. In the Apennine Peninsula, this juniper grows in the Central and Southern Apennines [10,11]. It occurs in the Alps [11,12] and, more rarely, in the Carpathians and in the Crimean Mountains [13,14]. Juniperus sabina also grows in Turkeyin the Manisa region and in the Black Sea region of the northern part of the country [15,16]. Within its range, J. sabina usually creates relatively small, scattered populations, often distant from each other.   Table 1. Symbols mean classes of the length of cone (CL) of populations. Light grey CL ≤ 6.00 mm. Dark grey 6.00 < CL ≤ 6.50 mm Black CL > 6.50 mm.

Biometry
The collected material was dried without pressing before measurement. Subsequently, 5 to 10 healthy and undamaged cones and 5 to 10 1-year-old lateral branchlets from each individual were measured. The traits for analyses were chosen based on previous examinations by different authors, performed for other junipers of the Sabina section [21,28,29]. To date, one of the features used to determine the size of the cones was cone diameter (CD), calculated as the average of two measurements of the width of the cone at its widest point, perpendicular to each other. We separately analysed the width (CW) and thickness of the cone (CTh), mainly because var. balkanensis is distinguished from var. sabina by having bi-lobed cones [24], indicating that the width of the cone should be greater than the thickness and that the CW/CTh ratio should be greater than 1. We also included the characteristic CD in the analyses; in total, we analysed 16 characteristics. Six of them were measured under the Delta Optical SZ-630 microscope: cone length (CL), cone width (CW), cone thickness (CTh), seed length (SL), seed width (SW) and thickness of the last ramification shoot with leaves (ST). The Delta Optical  Table 1. Symbols mean classes of the length of cone (CL) of populations. Light grey CL ≤ 6.00 mm. Dark grey 6.00 < CL ≤ 6.50 mm Black CL > 6.50 mm.
A sample from each population consisted of 10-30 individuals, depending on the size of the population, the presence of cone-bearing specimens and the availability of individual shrubs in a difficult to reach mountainous terrain, and therefore, in extremely small populations, fewer specimens were sampled. One sample from each individual comprised at least five properly developed, not injured, ripe cones and small parts of shoots. To avoid collecting the same genet, shrubs clearly delimited in the field, usually growing no closer than 30 m from each other, were treated as separate individuals.

Biometry
The collected material was dried without pressing before measurement. Subsequently, 5 to 10 healthy and undamaged cones and 5 to 10 1-year-old lateral branchlets from each individual were measured. The traits for analyses were chosen based on previous examinations by different authors, performed for other junipers of the Sabina section [21,28,29]. To date, one of the features used to determine the size of the cones was cone diameter (CD), calculated as the average of two measurements of the width of the cone at its widest point, perpendicular to each other. We separately analysed the width (CW) and thickness of the cone (CTh), mainly because var. balkanensis is distinguished from var. sabina by having bi-lobed cones [24], indicating that the width of the cone should be greater than the thickness and that the CW/CTh ratio should be greater than 1. We also included the characteristic CD in the analyses; in total, we analysed 16 characteristics. Six of them were measured under the Delta Optical SZ-630 microscope: cone length (CL), cone width (CW), cone thickness (CTh), seed length (SL), seed width (SW) and thickness of the last ramification shoot with leaves (ST). The Delta Optical DLTCamViewer software (Delta Optical, Nowe Osiny, Poland) was used to obtain digital images and to measure elements. Leaves were counted (LN) under the microscope on the last 5-mm long section of the shoot. Seed number (SN) was counted manually. Eight ratios were calculated with the measured/counted characteristics ( Table 2).

Statistical Analyses
Two types of data sets were investigated: (1) comprising all cases, i.e., each single measurement of the cone and shoot, and (2) consisting of the average values of individuals.
Descriptive statistics and properties were tested for non-transformed feature values, standardized values and logarithmically transformed values [30]. Unimodality and frequency distributions of values of features were checked using Shapiro-Wilk's test, and the distribution measures, skewness and kurtosis were calculated [31,32]. The descriptive statistics were arithmetic means, medians, modes, the lowest and highest values of characteristics, standard deviations, and coefficients of variation, which were compiled for all data, for data of each variety, each geographic region, and each population, as well as for sets of means of individuals of each of the above groups. Mutual correlations of characteristics were examined using Spearman's rank correlation coefficient as values of many characteristics were not normally distributed [33].
Further analyses were performed with the set of log-transformed values of individuals means. The homoscedasticity of variances of traits with normal distributions was verified with Brown-Forsythe's test [34]. Analysis of variance (ANOVA), followed by the Tukey's test, was performed for the traits with normal distributions and homoscedastic variations, and the Kruskal-Wallis test by ranks was done for the others [33]. These tests were conducted between populations, geographical regions and between two varieties and the groups of indefinite samples from the Carpathians and Tian Shan. The aim of these analyses was to delimit differences between the given groups, to find the best features to characterize them and to determine which groups differed by which characteristics. The results were illustrated via box-and-whiskers plots of chosen traits most important for differentiation, and shown in a common results table. The traits with normal distributions and homoscedastic variations were than used in the two-factor nested mixed ANOVA, for populations nested in regions.
The least mutually correlated, normally distributed characteristics, which, at the same time, most significantly influenced the sample differentiation in the ANOVA, were chosen for the discrimination analysis [35]. The analysis of discrimination from the multivariate statistics package of the STATISTICA 13 software (StatSoftPolska, Kraków, Poland) was performed to verify the level of differences and detect the geographic patterns, for the following groups: (1) two varieties: sabina and balkanensis and two indefinite groups: Carpathians and Tian Shan, (2) geographical regions, (3) all populations, (4) populations of each variety (two separate analyses) [32]. Results of these analyses were shown on the scatterplots between the first two discriminant variables [33]. Finally correlations between mean values of features and longitude, latitude and altitude of populations were calculated using the Spearman's rank correlation coefficient.
STATISTICA 13 software was used in all the above analyses.

Morphological Characteristics
The obtained values of features, in the set of all measurements, were not normally distributed, both for non-transformed and standardized data. Most of the distributions were right-skewed and leptokurtic, and the distributions of the measured characteristics were closer to normal compared to the ratios.  Table S1). Characteristics were moderately variable. The coefficient of variation (CV) of most characteristics took values between 10 and 20% (Table 2), with few exceptions. The ratio CW/CTh was the least variable, with the CV for all data equalling 4.90%, and in samples from 2.72% (SP6) to 8.17% (AU3). Ratio of leaves number/thickness of shoot (LN/ST) was the most variable feature, with a CV equalling 18.07% (from 8.25% in KYR1 to 22.79% in KYR2) (see Supplementary Material, Table S2).
Almost full correlations were found between the three features describing the width of the cone, namely CW, CTh and CD, and they were also highly correlated with CL, CL/CD, CD/SW and CD/ST, at p < 0.01. The parameter SL was significantly related to CL, whereas CW/CTh was not significantly correlated with any other parameter, and LN, ST and CL/SL were only correlated with one or two parameters each. Ratios were usually highly correlated to the characteristics involved; however CL/CD was not significantly correlated with CL, CL/SL was only weakly correlated with SL, and SL/SW was highly related to SW but weakly to SL.

Differences between Varieties, Regions and Populations
Comparisons between the varieties sabina and balkanensis and two indefinite groups of populations, namely Carpathians and Tian Shan, showed that all characteristics except seed width (SW) were responsible for differentiation, at p < 0.05. Most differences were found between the two varieties ( Table 3). The box-and-whisker graphs of means of varieties and means of populations showed the considerably greater variation of population means compared to the means of varieties, and the substantial variation of means of populations within regions (Figure 2a-c). According to the analyses performed to compare regions, the Iberian Peninsula and Crimea differed from other regions in terms of numerous features, whereas the Apennines and Tian Shan differed less ( Table 4). The comparisons between populations again showed that all characteristics differentiated them significantly, but the post hoc tests pointed out that most populations differed from others only in terms of a single feature, at p < 0.05. The greatest number of differences between pairs of samples was caused by CTh (177, (Table S3). Table 3. Characteristics differing significantly two varieties of J. sabina: sabina and balkanensis and two indefinite groups of populations: Carpathians and Tian Shan, according to results of Kruskal-Wallis' test; character codes as in Table 2.    Table 1.  Table 1.
Only the features: CTh, SL, ST, LN and CL/SL were at the same time weakly correlated, had the log-normal distributions and the strong impact on the differentiation of juniper groups. They were therefore used in the discrimination analyses. The grouping of individuals on the plot in the analysis performed for two varieties and two indefinite groups showed similar variations of the variety balkanensis and the Carpathian group, and the Tian Shan group covered a part of the variation of the variety sabina (Figure 3a). The differences between the two varieties were mostly influenced by ST and CTh. In the DA performed for regions, the grouping their means (Figure 3b) (Table 1). In square brackets−codes of characters the most strongly connected with the discriminant variable.  (Table 1). In square brackets-codes of characters the most strongly connected with the discriminant variable.
The analysis of correlations of feature values with geographic coordinates showed that only the number of leaves (LN) correlated positively with altitude (r = 0.57), and negatively with latitude (r = −0.67), with p < 0.01.   Table 1.

Morphological Characteristics
The results of our measurements were generally consistent with the data in monographs [4,17], although in a study by Amaral Franco [36], both cone width (4-6 mm) and shoot width (0.6-0.8 mm) were smaller than the values obtained in our research. These differences can be explained by the low number of specimens measured by the latter author. All characteristics selected in our study were useful, but the ratios showed several high correlations and were often non-parametric, limiting their use in multivariate analyses. Among the traits of cone width, namely CW, CTh and CD, the characteristic CTh played the greatest role in differentiating the analysed groups, but the differences among these three features were small.
The variability of our data was very similar to published data for J. excelsa [29] and J. phoenicea [23].

Intra-Species Differentiation
Our analyses showed statistically significant differences among populations and geographical regions, and between varieties of J. sabina, but the structure of this differentiation did not meet our expectations. The greatest differentiation was found among populations, and they did not cluster in any explicable way. Only populations of the two varieties analysed separately tended to show a weak geographical structure. This result, together with the statistically significant differences between the two varieties, indicates slightly different patterns of variation of these varieties. On the other hand, clustering of regions could imply the stronger effect of geographic location on differentiation, compared to the taxonomic status. Importantly, it should be noted that it was not possible to classify a given population as a variety based on the values of the features. The found differences between varieties can only be treated as a certain tendency, e.g., to a greater number of seeds for J. sabina.
The diagnostic value of the two features CW/CTh and SN, proposed by Adams et al. [24] as differentiating varieties, was not confirmed. The ratio of the two measures of cone width (CW/CTh) was the least variable characteristic with a biased distribution of its values, not significant in any comparisons. Seed number (SN) did not differ significantly between the two varieties, and its value was greater for the variety balkanensis (one to six seeds, with a median of three) than for the variety sabina (one to five seeds, with a median of two), contrary to the findings of Adams et al. [24]. As the comparisons made by Adams  Table 1.

Morphological Characteristics
The results of our measurements were generally consistent with the data in monographs [4,17], although in a study by Amaral Franco [36], both cone width (4-6 mm) and shoot width (0.6-0.8 mm) were smaller than the values obtained in our research. These differences can be explained by the low number of specimens measured by the latter author. All characteristics selected in our study were useful, but the ratios showed several high correlations and were often non-parametric, limiting their use in multivariate analyses. Among the traits of cone width, namely CW, CTh and CD, the characteristic CTh played the greatest role in differentiating the analysed groups, but the differences among these three features were small.
The variability of our data was very similar to published data for J. excelsa [29] and J. phoenicea [23].

Intra-Species Differentiation
Our analyses showed statistically significant differences among populations and geographical regions, and between varieties of J. sabina, but the structure of this differentiation did not meet our expectations. The greatest differentiation was found among populations, and they did not cluster in any explicable way. Only populations of the two varieties analysed separately tended to show a weak geographical structure. This result, together with the statistically significant differences between the two varieties, indicates slightly different patterns of variation of these varieties. On the other hand, clustering of regions could imply the stronger effect of geographic location on differentiation, compared to the taxonomic status. Importantly, it should be noted that it was not possible to classify a given population as a variety based on the values of the features. The found differences between varieties can only be treated as a certain tendency, e.g., to a greater number of seeds for J. sabina.
The diagnostic value of the two features CW/CTh and SN, proposed by Adams et al. [24] as differentiating varieties, was not confirmed. The ratio of the two measures of cone width (CW/CTh) was the least variable characteristic with a biased distribution of its values, not significant in any comparisons. Seed number (SN) did not differ significantly between the two varieties, and its value was greater for the variety balkanensis (one to six seeds, with a median of three) than for the variety sabina (one to five seeds, with a median of two), contrary to the findings of Adams et al. [24]. As the comparisons made by Adams et al. [24] were based on the first discoveries of the variety balkanensis in Greece and Bulgaria, their conclusion perhaps could only be applied to specimens of these two populations.
Although the differences between the varieties were weak, the indefinite samples from Tian Shan, according to our results, belonged to the variety sabina, whereas the variation of the Carpathian samples suggests their affinities to the variety balkanensis. These proposals, of course, require confirmation by genetic analyses, but the identification of samples from Tian Shan seems obvious. The balkanensis variety is the result of the ancient hybridisation between J. sabina and the ancestor of J. thurifera [24,27] and was found only on territories occupied by J. thurifera at present or in the past [28,37,38]. The identification of the Carpathian samples as var. balkanensis would be consistent with the proposal of Adams et al. [24].
Despite the differences discussed above, the geographical structure of the variation of J. sabina was weak and it could only be detected on the scatterplots of geographical regions and of populations of each variety. Besides, populations from the easternmost stands in Tian Shan, 4800 km from the Alps, did not show any specific features and grouped with the European populations. A similar geographic structure was obtained by Adams et al. [39] on the basis of research on essential oils.

Possible Reasons of the Weak Geographical Structure of J. sabina Morphology
The weak morphological differences between J. sabina var sabina and J. sabina var. balkanensis were partly expected, as the latter has been defined as morphologically cryptic [24]. Still, for several junipers in the Mediterranean region, a morphological intra-specific differentiation was recognised, which made us expect such variability in J. sabina. Intra-specific morphological differentiation, leading to discrimination between taxa, has been detected in the J. oxycedrus complex [40,41], J. thurifera [21] and J. phoenicea s.l. [23,42]. The variation of the Mediterranean junipers has been linked to their Tertiary and Pleistocene migrations [21,28,37,38] and long-term isolations [23,43]. The lack of geographic diversity of J. sabina, in view of the considerable distances between its populations, is puzzling. Based on the typical Mediterranean juniper species mentioned above, it is distinguished by its vast range, reaching Central Asia, and its mountainous nature [17]. Its dispersed type of the geographic range suggests the refugial characteristic of the distribution of the species, resembling those of alpine plant species, which were related to their migration after the end of the Pleistocene [44]. Morphological examinations of these Angiosperms based on vegetative traits revealed weak geographic signals, with climate as the possible cause of their intra-species variability [45][46][47]. However, a study of the complex Pinus mugo, based on vegetative features [48], as well as of P. sylvestris from the refugial areas, using both vegetative and generative morphological characteristics [49], found geographical differentiations. The homogeneity of the morphology of J. sabina could be the result of the long period of contact between nowadays distinct populations and their relatively late separation.
The species is extremely light-demanding, well adapted to a climate with dry and warm summers and relatively cold winters (with average temperatures of −6.5 to 32 • C), tolerant of snowfall, and resistant to drought [6,15,17,50,51]. Shrubs die when shaded by higher species of trees, such as spruce (personal observation made in the Alps and the Eastern Carpathians). It is therefore likely that the species occupied a larger territory during the Late Glacial and the Early Holocene. Because juniper species cannot be distinguished based on pollen and because of the scarcity of macrofossils, it is difficult to trace the history of these species [52][53][54]. The widespread occurrence of the species could have been promoted by the dryness of the climate during the Younger Dryas and again in the Preboreal and Boreal periods [52], with an increase in seasonality and with the highest precipitation in winter [55], and by the significantly lower treeline in the mountains [53,56]. The pressure of deciduous forests developing during the Atlantic period could have limited the occurrence of junipers, which was found, for example, on the Iberian Peninsula [52].
The survival of the species was possible because of the plasticity of its response to changing climatic conditions [57]. The subsequent drying of the climate once again favoured the spread of juniper, although human activity also increased during this period, especially pastoralism, which could have had an impact on juniper [55].
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/d13100470/s1, Table S1: Normality statistics of characteristics (characteristic codes- Table 2) for data, means of individuals (MI) and logarithms (log), of varieties and regions (Table 1), Table  S2: Mean values (M) and coefficient of variations (CV) of characteristics (codes as in Table 2) for all populations (acronyms-Table 1), Table S3: Results of two-factor nested mixed analysis of variance between populations and regions, for characters meeting assumptions.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available from the corresponding author upon request.