Soil Requirements, Genetic Diversity and Population History of the Juniperus sabina L. Varieties in Europe and Asia

: Trees and shrubs belonging to the genus Juniperus L. are pivotal species in arid and semiarid ecosystems in the Northern Hemisphere. However, unfavourable phenomena are observed in their populations due to global warming. We aimed to investigate the soil requirements, genetic diversity and population history of Juniperus sabina L. from Europe, Georgia, and Kyrgyzstan. Genetic resources were evaluated in 16 populations using nuclear microsatellites, while past demographic events were described based on the chloroplast DNA haplotypes. Seven chemical parameters in 36 soil samples from the European range of J. sabina were compared. In the studied area, three distinct phylogenetic lineages corresponding to different varieties of J. sabina , namely var. sabina , var. balkanensis , and the Asian variety, were revealed. Unimodal mismatch distributions and significantly negative Tajima’s D and Fu’s Fs parameters indicated that the sabina and balkanensis varieties underwent a population expansion. Microsatellite variation was moderate, potentially influenced by inbreeding, clonal propagation, and limited gene flow between populations. Bayesian clustering revealed five genetic groups. Compared to var. sabina , the balkanensis variety occupies areas with significantly higher potassium content in the soil, which probably mitigates the adverse effects of drought in its localities.


Introduction
Rising global temperatures, increased water scarcity, and extreme weather events threaten Earth's biodiversity.Population survival in warming and drying conditions is related to the fecundity of individuals, and adaptive capacity, defined by the level of genetic diversity, phenotypic plasticity, and connectivity with other populations [1,2].Among the plant taxa affected by climate deterioration are trees and shrubs of the genus Juniperus L. that inhabit arid and semiarid ecosystems in the Northern Hemisphere [3][4][5].In the Przewalski's juniper (Juniperus przewalskii Kom.) populations in the northeastern Tibetan Plateau, annual ring narrowing and missing rings in the tree trunks were observed due to four extreme droughts in the twentieth century [6].Up to 85% of medium-to large-sized Ashe juniper (Juniperus ashei J. Buchholz) trees may not survive a severe drought in Central Texas [7].
The history of the juniper dates back to the beginning of the Cenozoic [5].Its diversification probably started in the Eocene and proceeded at different rates throughout the Tertiary.This process was usually linked to global and local climate changes [4,5], involving hybridisation and allo-polyploidisation events [8][9][10].The genus is currently divided into three monophyletic sections: Caryocedrus, Juniperus and Sabina [3].Savin juniper (Juniperus sabina L.) is a Eurasian representative of the Sabina section (Figure 1).In Europe, the species Forests 2024, 15, 866 2 of 18 includes two varieties: diploid var.sabina and tetraploid var.balkanensis R.P. Adams and A.N. Tashev [3].The tetraploid variety has arisen via crossing between the maternal var.sabina and an ancestor of paternal tetraploid Spanish juniper (Juniperus thurifera L.) [3,8,11].Four hybridisation scenarios were proposed [8].One proposal involves fertilisation of unreduced female gametes by haploid pollen, resulting directly in tetraploid progeny, and the other three possibilities presume the existence of a 'triploid bridge' before the tetraploid genome was formed.The 'triploid bridge' hypotheses assume that triploid interspecific hybrids produced reduced gametes.According to all scenarios, backcrossing with maternal plants was needed [8].
Forests 2024, 15, 866 2 of 17 into three monophyletic sections: Caryocedrus, Juniperus and Sabina [3].Savin juniper (Juniperus sabina L.) is a Eurasian representative of the Sabina section (Figure 1).In Europe, the species includes two varieties: diploid var.sabina and tetraploid var.balkanensis R.P. Adams and A.N. Tashev [3].The tetraploid variety has arisen via crossing between the maternal var.sabina and an ancestor of paternal tetraploid Spanish juniper (Juniperus thurifera L.) [3,8,11].Four hybridisation scenarios were proposed [8].One proposal involves fertilisation of unreduced female gametes by haploid pollen, resulting directly in tetraploid progeny, and the other three possibilities presume the existence of a 'triploid bridge' before the tetraploid genome was formed.The 'triploid bridge' hypotheses assume that triploid interspecific hybrids produced reduced gametes.According to all scenarios, backcrossing with maternal plants was needed [8].The range of var.sabina extends from Spain in the west to Georgia in the east, whereas var.balkanensis is found in the Iberian, Apennine and Balkan Peninsulas and western Turkey [11][12][13][14][15].The varieties sabina and balkanensis co-occur only in one locality in Spain, but hybrids were not found there [9].The almost non-overlapping distribution of the two juniper varieties results from their different climatic requirements [14].The sabina variety inhabits areas with substantial precipitation and the lowest average annual temperatures, while var.balkanensis is adapted to the highest annual temperatures and intermediate precipitation levels compared to parental taxa.Plant survival in territories with water scarcity requires the evolution of different adaptive traits, e.g., early flowering, rapid development, transpiration reduction, trichomes formation or changes in the root system [16], and references therein.In the semiarid parts of China, the J. sabina specimens form vertical root systems homogenously distributed from the soil surface to groundwater levels [17].This feature enables moisture extraction from different layers within the soil profile.It is commonly known that the adverse effects of drought stress in plant crops can be mitigated by the potassium (K) supply [16,18].Potassium is fundamental to plant life as it serves as an essential macronutrient in numerous physiological processes, including protein synthesis, enzyme activation, phloem transport, photosynthesis, stomatal regulation, and osmoregulation [18][19][20].It was shown that the midday leaf temperature of potatoes (Solanum tuberosum L.) growing under drought stress was significantly lowered due to increased leaf water potential from the K supply [21].The balkanensis variety of J. sabina exhibits a substantial drought tolerance [15]; however, factors enhancing its ability to cope with the water deficit have not yet been recognised.The range of var.sabina extends from Spain in the west to Georgia in the east, whereas var.balkanensis is found in the Iberian, Apennine and Balkan Peninsulas and western Turkey [11][12][13][14][15].The varieties sabina and balkanensis co-occur only in one locality in Spain, but hybrids were not found there [9].The almost non-overlapping distribution of the two juniper varieties results from their different climatic requirements [14].The sabina variety inhabits areas with substantial precipitation and the lowest average annual temperatures, while var.balkanensis is adapted to the highest annual temperatures and intermediate precipitation levels compared to parental taxa.Plant survival in territories with water scarcity requires the evolution of different adaptive traits, e.g., early flowering, rapid development, transpiration reduction, trichomes formation or changes in the root system [16], and references therein.In the semiarid parts of China, the J. sabina specimens form vertical root systems homogenously distributed from the soil surface to groundwater levels [17].This feature enables moisture extraction from different layers within the soil profile.It is commonly known that the adverse effects of drought stress in plant crops can be mitigated by the potassium (K) supply [16,18].Potassium is fundamental to plant life as it serves as an essential macronutrient in numerous physiological processes, including protein synthesis, enzyme activation, phloem transport, photosynthesis, stomatal regulation, and osmoregulation [18][19][20].It was shown that the midday leaf temperature of potatoes (Solanum tuberosum L.) growing under drought stress was significantly lowered due to increased leaf water potential from the K supply [21].The balkanensis variety of J. sabina exhibits a substantial drought tolerance [15]; however, factors enhancing its ability to cope with the water deficit have not yet been recognised.The sabina and balkanensis varieties are morphologically indistinguishable [22], but they are genetically distinct respecting chloroplast (cp) DNA and both SilicoDArT and SNP loci obtained by whole-genome sequencing [9,11,14].The SilicoDArT-and SNPbased analyses revealed relatively low levels of genetic diversity of both varieties, but this investigation involved a small number of individuals per population [14].In var.balkanensis populations from the Balkan Peninsula, only three individuals per population were studied using the inter-simple sequence repeat (ISSR) markers [15].In the present study, the nuclear microsatellite and cpDNA markers were used to analyse genetic diversity and differentiation and the population history of the savin juniper in Europe and two Asian countries, Georgia and Kyrgyzstan.As the J. sabina range is strongly fragmented and the species populates ecosystems that are impacted by climate deterioration, we specifically aimed to: (1) determine genetic variation in the populations, (2) infer potential population structuring, (3) examine a possible relationship between the soil chemical parameters and the distribution of sabina and balkanensis varieties, and (4) describe demographic changes in the populations studied.
Forests 2024, 15, 866 3 of 17 The sabina and balkanensis varieties are morphologically indistinguishable [22], but they are genetically distinct respecting chloroplast (cp) DNA and both SilicoDArT and SNP loci obtained by whole-genome sequencing [9,11,14].The SilicoDArT-and SNPbased analyses revealed relatively low levels of genetic diversity of both varieties, but this investigation involved a small number of individuals per population [14].In var.balkanensis populations from the Balkan Peninsula, only three individuals per population were studied using the inter-simple sequence repeat (ISSR) markers [15].In the present study, the nuclear microsatellite and cpDNA markers were used to analyse genetic diversity and differentiation and the population history of the savin juniper in Europe and two Asian countries, Georgia and Kyrgyzstan.As the J. sabina range is strongly fragmented and the species populates ecosystems that are impacted by climate deterioration, we specifically aimed to: (1) determine genetic variation in the populations, (2) infer potential population structuring, (3) examine a possible relationship between the soil chemical parameters and the distribution of sabina and balkanensis varieties, and (4) describe demographic changes in the populations studied.
Table 1.Genetic characteristics of J. sabina populations.N PL -ploidy level, N-number of samples, N G -number of distinct genetic individuals, A-mean number of alleles per locus, H O -observed heterozygosity, H S -gene diversity, G IS -inbreeding coefficient, P(HWE)-p value for HWE, *-statistically significant values after sequential Bonferroni's correction.The population acronyms are the same as in [14].

Testing and Choice of the Nuclear Microsatellite Loci
In total, 15 nuclear microsatellite loci were tested on the set of 32 individuals from different populations.These were six primer pairs developed for the J. sabina loci (Sabv6, Sabv10, Sabv13, Sabv15, Sabv16, Sabv17) [23], four pairs designed for the common juniper (Juniperus communis L.) loci (Jc031, Jc032, Jc037, Jc166) [24,25], and five primer pairs designed for the Phoenicean juniper (Juniperus phoenicea ssp.turbinate) loci (Junpho_005458, Junpho_006365, Junpho_064797, Junpho_084596, Junpho_127300) [26].However, the loci Jc031, Jc037, and Junpho_064797 failed to amplify in all samples; seven loci (Sabv10, Sabv13, Sabv17, Jc032, Junpho_006365, Junpho_084596, Junpho_127300) showed amplification in some samples only, and loci Sabv16 and Junpho_005458 revealed multi-banding patterns in the tetraploids.For the above reasons, 12 loci were discarded at this pilot step.The primers for the Jc166 locus did not succeed in generating amplification products in the samples from the KYR population, but we decided to include this locus in the study because of the good quality of the products in all the other samples.Finally, Sabv6 and Sabv15 loci were amplified in the multiplex, while Jc166 locus was amplified independently.Please refer to Table S1 for the sequences of the nuclear microsatellite primers tested in this study.
All PCR conditions consisted of an initial hot-start denaturation at 95 • C for 15 min, followed by 38 for Sabv6 and Sabv15 or 34 for Jc166 cycles at 94 • C for 30 s, 55 • C for 45 s, 72 • C for 50 s and a final extension step at 72 • C for 7 min.PCRs were performed in a volume of 5 µL containing 2 µL of diluted template DNA, 1.68 µL of StartWarm HS-PCR Mix (A&A Biotechnology, Gda ńsk, Poland), 0.99 µL of RNase-free water (A&A Biotechnology), and 0.33 µL of mixed primers (0.2 µM of each primer).PCRs were performed using a Labcycler Basic (SensoQuest, Göttingen, Germany).The microsatellite fragments were separated in an ABI PRISM 3100 Genetic Analyser (Applied Biosystems, Foster City, CA, USA) with a Gene Scan-500 LIZ size standard (Applied Biosystems).Genotypes were scored with the help of GeneMapper 4.0 software (Applied Biosystems).The initial screening of genotypes involved individuals from the confirmed diploid populations (SP7, AU3, AU4, RO3, IT1, PL, GEO) only.Then, the allele copy number ('dosage') was established in the remaining individuals based on peak height ( [27], and references therein) made independently by two persons.

Soil Data
The soil data were obtained from the European Soil Data Centre (ESDAC) (https: //esdac.jrc.ec.europa.eu/;accessed on 16 February 2024).The dataset was collected from the Land Use/Cover Area Frame Statistical Survey (LUCAS) 2018 [28,29].The data file contained chemical properties for the soil samples collected around Europe.These were: total nitrogen content (N; g × kg −1 ), extractable phosphorus (P; mg × kg −1 ) and potassium (K; mg × kg −1 ) contents, carbonates content at depth 0-20 cm (CaCO 3 ; g × kg −1 ), the electrical conductivity (EC; mS × m −1 ), organic carbon content at depth 0-20 cm (OC; g × kg −1 ), and pH measured in a suspension of soil in water.Out of 18,984 records, we selected 36 data points most geographically adjacent to the J. sabina populations under study.These were 20 points for var.sabina and 16 points for var.balkanensis.

Statistical Analyses 2.4.1. Genetic Diversity of the Nuclear Microsatellites
The software package GenoDive Version 3.06 [30], comprising a correction for unknown dosage of alleles in the nuclear microsatellite loci in multiple analyses, was used to assess the genetic diversity and differentiation parameters within and between J. sabina populations.The number of distinct genotypes (genets) was established by deciding that missing data would not be counted.After the removal of duplicate genotypes (clones; Table 1), the presence of null alleles, stutter peaks and allele dropout was checked in the diploid populations (SP7, AU3, AU4, RO3, IT1, PL, GEO, KYR) using the MICRO-CHECKER 2.2.3 [31].The frequency of null alleles and p-values for heterozygote deficiency were calculated using 10,000 Monte Carlo randomisations implemented in the ML-NullFreq software (version 2006) [32].The number of alleles (N A ), observed heterozygosity (H O ), gene diversity (H S ), inbreeding coefficient (G IS ) and the deviations from Hardy-Weinberg equilibrium (HWE) were assessed for each locus in GenoDive.The mean number of alleles per locus (A), H O , H S , G IS , and the departures from HWE were calculated at the population level.Significant departures from HWE at both locus and population levels were tested using 10,000 permutations.

Genetic Differentiation and Population Structure
The ρ (Rho) statistic, which is comparable between the different ploidy levels and is independent of the double reduction rate, implemented in GenoDive was used to estimate genetic differentiation between pairs of populations using the microsatellites [33].The significance of ρ was tested by running 10,000 permutations, and the results obtained were corrected using the sequential Bonferroni correction [34].STRUCTURE 2.3.4 software [35] was used to discover the most likely number of independent genetic clusters (K) performing ten independent runs for K values ranging from 1 to 16 with 500,000 iterations and 50,000 burn-in periods.All J. sabina genotypes in the input file were coded as tetraploid; however, two lacking chromosome sets were assigned as 'missing data' for diploid individuals.The Bayesian clustering investigation was based on two models: LOCPRIOR (locality information) and admixture with correlated allele frequencies.Due to uneven sample sizes, Puechmaille's approach [36] was chosen and four statistical measures, MedMed (median of medians), MedMean (median of means), MaxMed (maximum of medians), and MaxMean (maximum of means), were calculated to indicate the optimal K value using StructureSelector software (version 2018) [37].The population was assigned to a cluster if the membership coefficient (=ancestry) to that cluster was higher than 0.5 (threshold value) [36].The STRUCTURE results were summarised and visualised using the CLUMPAK software (version 2018) [38].

CpDNA Marker Analyses
All cpDNA analyses were based on the haplotypes described previously for 108 J. sabina individuals collected in 17 populations from Europe and Asia [14].Haplotypes were established based on the four intergenic spacers (trnL-trnF, trnD-trnT, trnS-trnG, petN-psbM), and the total length of combined cpDNA sequence was 3024 bp (for the GenBank accession numbers and other details see [14]).The minimum spanning network (MSN) was constructed using a statistical parsimony network approach implemented in TCS 1.21 [39] to establish the relationships among the cpDNA haplotypes.In this reconstruction, each indel was treated as a fifth state.
To infer the past demographic events in the J. sabina populations, the empirical mismatch distribution (distribution of nucleotide-site differences between pairs of the sequences) was compared with theoretical distribution under the population growth-decline model implemented in DnaSP version 6 [40].Tajima's D [41] and Fu's Fs test statistic [42] were calculated assuming the neutral evolution model to check access of rare mutations as evidence of population expansion.The statistical significance of D and F parameters were tested using the coalescent algorithm with 10,000 bootstrap replicates.All the above calculations were conducted in the total sample (108 sequences) and three subsamples: all diploids (60), European diploids and GEO (45), and tetraploids (48).

Analyses of the Soil Parameters
Two-sample randomisation tests (test statistic-median difference; 10,000 randomisations) implemented in RundomPro 3.14 [43] were conducted to compare seven soil parameters between the localities of var.sabina and var.balkanensis.Next, the p values were corrected using the sequential Bonferroni correction [34].To test the potential impact of soil parameters (explanatory variables) on the distribution of var.sabina and var.balkanensis (response variables) in Europe, redundancy analysis (RDA) with 9999 permutations using PAST version 4.11 software [44] was carried out.Before this analysis, the values of soil parameters were standardised using the Z score method.

Diversity of Microsatellite Loci and Populations
The highest number of alleles (16) was found in the Jc166 locus, with a mean of 13.3 per locus (Table 2).In turn, the highest observed heterozygosity (H O = 0.73) and gene diversity (H S = 0.75) were noted in Sabv6, with mean values of 0.44 and 0.59 for H O and H S , respectively.The inbreeding coefficient (G IS ) ranged from 0.02 in Sabv6 to 0.44 in Jc166.Statistically significant deviations from HWE were revealed in the loci Sabv15 and Jc166.Among 335 genotyped samples, 264 genets of J. sabina were identified (Table 1).The KYR specimens were genotyped using two loci, as the primers for the Jc166 locus failed to generate amplification products in those individuals.Among possible genotyping errors, the null alleles in locus Jc166 showed statistically significant frequencies ranging from 0.19 to 0.41 were detected in the SP7, AU4, RO3, and GEO populations.Null alleles were also revealed in the locus Sabv15 in IT1 (0.18, p = 0.0081) and RO3 (0.16, p = 0.0001) populations.The number of alleles per locus ranged from 2.3 in PL to 6.0 in GR, H O was between 0.34 (KYR) and 0.61 (IT3), while H S ranged from 0.31 in PL to 0.72 in IT3 and IT4.Except for the PL population, where no inbreeding was revealed (G IS = −0.22), the remaining G IS values were high, and they ranged from 0.14 in BG to 0.38 in RO3, with a mean of 0.25.The significant departure from HWE was noted in 11 populations, excluding the AU3, AU4, IT1, PL, and KYR.The total sample also significantly deviated from HWE.

Genetic Differentiation and Population Structure
Out of 120 pairwise comparisons, statistically significant genetic differentiation (ρ) values were revealed between 59 population pairs (Table 3).The highest differentiation was between SP7 and PL (ρ = 0.408, p < 0.0001); the population pairs AU4 and NM, RO3 and IT1, RO3 and IT3, RO3 and CRO, IT1 and KYR, and IT3 and CRO were not genetically differentiated.Puechmaille's statistical measures [36] indicated five genetic clusters in the area studied (Figure S1); however, nine populations only achieved the mean membership coefficient ≥ 0.5, which allowed us to include them univocally in one of the groups.The first group (orange colour) included SP7 and KYR samples (Figure 3).The second cluster (blue) consisted of five populations: AU3, AU4, NM, PL, and GEO.The third (dark red) and fourth (green) clusters involved single populations, RO1 and BG, respectively.The fifth cluster was spurious, i.e., a cluster involving populations with a mean ancestry coefficient lower than a threshold value.) inferred by Puechmaille's approach [36] with nuclear microsatellite markers' dataset.Each colour corresponds to one of five genetic groups.

CpDNA-Based Analyses of Population History
The minimum spanning network (MSN) involving 19 cpDNA haplotypes (H1-H19), described previously [14], revealed three haplogroups deeply differentiated from each other (Figure 4).The first haplogroup comprised a widespread H3 haplotype and 11 additional ones discovered in the sabina samples from Europe and the GEO locality from Asia.The second haplogroup included the common H14 and four related haplotypes from the var.balkanensis populations.The third haplogroup consisted of H18 and H19 from the KYR sample from Asia.The haplotypes H3 and H14 differed by 44 mutation steps from one another.There were ten mutations between European H3 and Kyrgyz H18 haplotypes and 35 between H18 and H14.[36] with nuclear microsatellite markers' dataset.Each colour corresponds to one of five genetic groups.

CpDNA-Based Analyses of Population History
The minimum spanning network (MSN) involving 19 cpDNA haplotypes (H1-H19), described previously [14], revealed three haplogroups deeply differentiated from each other (Figure 4).The first haplogroup comprised a widespread H3 haplotype and 11 additional ones discovered in the sabina samples from Europe and the GEO locality from Asia.The second haplogroup included the common H14 and four related haplotypes from the var.balkanensis populations.The third haplogroup consisted of H18 and H19 from the KYR sample from Asia.The haplotypes H3 and H14 differed by 44 mutation steps from one another.There were ten mutations between European H3 and Kyrgyz H18 haplotypes and 35 between H18 and H14.ditional ones discovered in the sabina samples from Europe and the GEO locality from Asia.The second haplogroup included the common H14 and four related haplotypes from the var.balkanensis populations.The third haplogroup consisted of H18 and H19 from the KYR sample from Asia.The haplotypes H3 and H14 differed by 44 mutation steps from one another.There were ten mutations between European H3 and Kyrgyz H18 haplotypes and 35 between H18 and H14.In the total sample, positive and statistically non-significant Tajima's D (2.680) and Fu's Fs (17.481) suggested a deficiency of rare mutations over what would be expected under the neutral model (Table 4).The values of Tajima's D and Fu's Fs statistics were  [14].CR-Crimea; the remaining population names are according to Table 1.Blue colour indicates haplotypes of var.sabina, red-var.balkanensis, and green-the Asian variety from Kyrgyzstan.
In the total sample, positive and statistically non-significant Tajima's D (2.680) and Fu's Fs (17.481) suggested a deficiency of rare mutations over what would be expected under the neutral model (Table 4).The values of Tajima's D and Fu's Fs statistics were negative in three subsets of data; however, they were statistically significant in the group comprising diploids from Europe and GEO (D = −2.186,p = 0.000; Fs = −3.516,p = 0.004) and in the set of tetraploids (D = −1.107,p = 0.004; Fs = −1.602,p = 0.006).Negative values of both parameters indicate an excess of rare mutations relative to expectation.The mismatch distributions for the total dataset were multimodal (Figure 5A); this was bimodal for the set containing all diploid populations (Figure 5B).European diploids with the GEO sample and tetraploid populations revealed skewed unimodal mismatch distributions (Figure 5C,D).

Soil Parameters
Substantially higher content of K (median difference-101.2,p = 0.002) was foun the soils overgrown by var.balkanensis compared to the sabina lands (Figure 6).The va of remaining soil parameters were not significantly different between the two grou junipers.The initial RDA revealed strong correlations between different soil fa

Soil Parameters
Substantially higher content of K (median difference-101.2,p = 0.002) was found in the soils overgrown by var.balkanensis compared to the sabina lands (Figure 6).The values of remaining soil parameters were not significantly different between the two groups of junipers.The initial RDA revealed strong correlations between different soil factors (overlapping the arrowheads of particular explanatory variables); thus, four parameters (P, K, pH, CaCO 3 ) were used in the final analysis.The RDA revealed that K concentrations had the most critical influence on the distributions of savin juniper varieties, as the arrowhead of that explanatory variable was the longest (Figure 7).K content was positively correlated with var.balkanensis and negatively with var.sabina, as the angle between the K arrowhead and former and latter juniper varieties was acute and obtuse, respectively.Additionally, var.sabina locations were positively related to P content, while higher pH was found in the var.balkanensis locations.Consequently, the sabina and balkanensis data points were somewhat separated in the RDA plot (F = 6.818, p = 0.001).The first eigenvalue (0.23772) implied that the total variance revealed in RDA (46.8%) was explained by the first axis because the second eigenvalue was almost zero (9.3435 × 10 −33 ).
ditionally, var.sabina locations were positively related to P content, while higher p found in the var.balkanensis locations.Consequently, the sabina and balkanensis data were somewhat separated in the RDA plot (F = 6.818, p = 0.001).The first eige (0.23772) implied that the total variance revealed in RDA (46.8%) was explained first axis because the second eigenvalue was almost zero (9.3435 × 10 −33 ).

Chemical Soil Differences in the European Range of J. sabina
The salient result of our study was that the occurrence of the sabina and balkanensis varieties in Europe was impacted not only by the climate conditions [14] but also by the different chemical properties of the soils (RDA; F = 6.818, p = 0.001).The var. balkanensis overgrows areas with significantly higher extractable K content than var.sabina (median difference-101.2,p = 0. 0.002).Potassium regulates many physiological functions of plant cells and organs, e.g., synthesis of proteins, enzyme activation, phloem transport, photosynthesis, stomatal regulation, and osmoregulation; thus, it is one of the critical macroelements in plant life [18][19][20].Evidence indicates that this nutrient supports the survival of plants under different biotic and abiotic stress conditions [16,20].Plant crops with a proper supply of K were shown to be less susceptible to insect infestation and microbial infections than K-deficient cultivations [20].Applying K mitigates the detrimental impacts of water deficit, enabling plants subjected to drought stress to thrive and produce fruit [16,18].Adequate K resources in soils assist in sustaining plant resistance to drought through diverse mechanisms: increase in root length and root density, stability and regulation of the osmotic potential and hydraulic conductivity of cell membranes, formation of osmotic adjustment ability, regulation of turgor within the guard cells during stomatal movement, and the reduction of reactive oxygen species (ROS) production [16,18,20].It was revealed that the rapeseed (Brassica napus L.) K-supplied cultivars were stimulated to release a more significant amount of organic acids in the root, which enhanced nutrient acquisition and utilisation under drought conditions [18].
The most important factor responsible for the K variations in the soils is climate, but topographical conditions, clay mineralogical composition, soil properties, and land use also influence the K levels [45,46].In Europe, the highest K content was detected in the soils characterised by the highest contribution of clay [47].Clay soils predominate in the southern part of the continent, in the Mediterranean Basin [45].Low K concentrations are found in the sandy soils in northeastern and northern European countries [45,47].Compared to the var.sabina localities, the var.balkanensis populations occur under warmer and drier climate conditions [14]; thus, the substantially higher K concentrations in its  The salient result of our study was that the occurrence of the sabina and balkanensis varieties in Europe was impacted not only by the climate conditions [14] but also by the different chemical properties of the soils (RDA; F = 6.818, p = 0.001).The var. balkanensis overgrows areas with significantly higher extractable K content than var.sabina (median difference-101.2,p = 0. 0.002).Potassium regulates many physiological functions of plant cells and organs, e.g., synthesis of proteins, enzyme activation, phloem transport, photosynthesis, stomatal regulation, and osmoregulation; thus, it is one of the critical macroelements in plant life [18][19][20].Evidence indicates that this nutrient supports the survival of plants under different biotic and abiotic stress conditions [16,20].Plant crops with a proper supply of K were shown to be less susceptible to insect infestation and microbial infections than K-deficient cultivations [20].Applying K mitigates the detrimental impacts of water deficit, enabling plants subjected to drought stress to thrive and produce fruit [16,18].Adequate K resources in soils assist in sustaining plant resistance to drought through diverse mechanisms: increase in root length and root density, stability and regulation of the osmotic potential and hydraulic conductivity of cell membranes, formation of osmotic adjustment ability, regulation of turgor within the guard cells during stomatal movement, and the reduction of reactive oxygen species (ROS) production [16,18,20].It was revealed that the rapeseed (Brassica napus L.) K-supplied cultivars were stimulated to release a more significant amount of organic acids in the root, which enhanced nutrient acquisition and utilisation under drought conditions [18].
The most important factor responsible for the K variations in the soils is climate, but topographical conditions, clay mineralogical composition, soil properties, and land use also influence the K levels [45,46].In Europe, the highest K content was detected in the soils characterised by the highest contribution of clay [47].Clay soils predominate in the southern part of the continent, in the Mediterranean Basin [45].Low K concentrations are found in the sandy soils in northeastern and northern European countries [45,47].Compared to the var.sabina localities, the var.balkanensis populations occur under warmer and drier climate conditions [14]; thus, the substantially higher K concentrations in its localities than in the sites of var.sabina can reduce the deleterious effect of drought stress on this juniper.It can be expected that the balkanensis shrubs, with access to adequate K content, can develop long and dense root systems, enabling them to capture soil water more efficiently and from deeper levels in the soil profile than the sabina individuals.Such a root system was described in the J. sabina individuals inhabiting the Mu Us Sandy Land in China [17].
RDA in our study also showed a tendency for var.balkanensis to appear in soils with a higher pH and lower P content than var.sabina; however, these results were not statistically significant.We suppose that higher median value of pH in the var.balkanensis populations (6.31) compared to the var.sabina sites (5.61) could be related to higher K-contents in the Mediterranean Basin [47].In the experiments involving various rapeseed cultivars, Xu et al. demonstrated that K-fertilizer had a positive impact on both soil pH and available phosphorus levels, consequently enhancing the diversity of the microbial community [18].The improved soil environment could assist plants in maintaining resistance against water stress [18].The extractable P resources in all J. sabina locations in Europe seem to be low (up to 70 mg × kg −1 ), although phosphorus is an essential macro-element for plant growth and reproduction [19].It was revealed that bioavailable P content was negatively correlated with soil pH in the alpine meadows in Tian Shan Mts.[48], consistent with our findings.The primary morphological adaptations that enable plants to efficiently utilise low P levels include the elongation of root length and the proliferation of lateral roots [49].
Possible differences in adaptations of the J. sabina varieties to the environmental conditions have not been analysed so far.Morphological investigation indicated that the cone thickness and shoot thickness of var.balkanensis were greater than those found in var.sabina, but the classification of a population into one or another variety failed based on the individual features [22].It would be worthwhile to study how var.balkanensis is morphologically and physiologically adapted to low precipitation and high temperatures in order to assess its potential survival under climate change.

Genetic Diversity of the J. sabina Populations
The nuclear microsatellite variation in 16 J.sabina populations spread over the area from Spain to Kyrgyzstan was moderate, with mean A = 4.6, H O = 0.44, and H S = 0.59.However, these values were higher than those estimated in 11 J. sabina populations in China (A = 3.64, H O = 0.402, H E = 0.456) [50].The most genetically diverse populations in our study were those of var.balkanensis: IT3 (A = 5.6, H O = 0.61, H S = 0.72) and GR (A = 6.0,H O = 0.60, H S = 0.71), while the lowest values of genetic variation parameters were revealed in the smallest PL stand of var.sabina (A = 2.3, H O = 0.38, H S = 0.31).Several factors can explain the limited genetic diversity observed in the savin juniper populations.First, generative reproduction is likely severely limited in many localities.Only 2.5% of the J. sabina seeds in southern Mongolia were viable [51].In the Polish population of var.sabina in the Pieniny Mts., the abundance of filled seeds was rare [52].Second, the savin juniper populations are dominated by clonal plants.Re-rooting clonal fragments, a primary mode of reproduction in the poor-soil desert steppes and mountains [53], was a likely reason for lowered genetic variation of J. sabina in northern China [54].Our study revealed the highest number of clones in the PL (68%) and KYR (54%) populations.The PL accessions were cultivated from seeds in the Kórnik Arboretum (Institute of Dendrology, Polish Academy of Sciences) [52]; thus, at least some of them can share a common maternal origin.In turn, two microsatellite loci were amplified only in the KYR population, likely falsifying the actual number of genetically identical specimens.
The third reason for the restricted genetic resources of the J. sabina populations is inbreeding.We found the positive values of the inbreeding coefficient in 15 localities (G IS = 0.14-0.38),while significant deviations from HWE were noted in 11 populations.The SNP markers also suggested a substantial level of inbreeding in the samples of 94 savin juniper individuals (F IS = 0.36) [14].Based on the nuclear microsatellites, crossing between close relatives was detected in 11 Chinese populations of J. sabina (F IS = 0.010-0.282)[50].It was implied that the substantial contributions of null alleles (0%-0.335%) in two Chinese savin juniper populations could result from inbreeding, as positive F IS values (0.050-0.699) were revealed in seven out of eight loci [23].We found the highest frequencies of null alleles in the Jc166 locus, but such alleles were also noted in the Sabv15 locus in the IT1 and RO1 populations.The null alleles were also detected in the J. communis populations in Germany (12%-41%) [25].It is not excluded that the presence of null alleles in different Juniperus studies could be simply a consequence of the cross-species/variety amplification problems of the nuclear microsatellites.Among 31 microsatellite loci developed for J. communis and J. przewalskii tested in the Greek juniper (Juniperus excelsa M. Bieb.) populations, the amplification products were obtained in seven loci only, but four loci were monomorphic [24].
Numerous statistically significant results obtained from the pairwise comparisons of genetic differentiation indicate that increased genetic diversity of J. sabina populations through gene exchange is unrealistic.Limitation in gene flow in the studied area was expected because the populations inhabit different mountain ranges, resulting in severe fragmentation of the savin juniper's range.Analyses based on the SilicoDArT and SNP markers also indicated the limitation in gene flow among the J. sabina localities and significant isolation by distance pattern [14].Based on the statistically significant Rho values in the pairwise comparisons of genetic differentiation in the microsatellite analysis, SP7 and PL populations of var.sabina seemed to be most different from the remaining populations.The SP7 locality is approximately 1030 km away from the nearest IT1 population, preventing gene flow between them.The PL population from the northern edge of the species range can be the smallest among the studied populations; consequently, it can be subjected to strong genetic drift, which drives this locality to differentiation.
The STRUCTURE analysis using the microsatellites revealed a clustering pattern (K = 5) that was challenging to interpret.On the one hand, the KYR population was grouped with SP7 from Spain.On the other hand, NM representing var.balkanensis was included in the cluster with var.sabina populations AU3, AU4, PL, and GEO.Among the var.balkanensis localities, the RO1 from the Carpathians and BG from the Balkans were clearly distinct from each other and from the remaining populations.In the STRUCTURE based on SNPs, KYR and NM also clustered with other var.sabina populations [14].While the challenges in delineating population structure within the J. sabina range using microsatellites may partly stem from the limited number of loci, we believe that a more significant factor is the substantial contribution of var.sabina nuclear genes in the var.balkanensis populations [11,14].Morphological analyses showed a striking similarity between the KYR individuals and var.sabina from Europe [22].The morphological differences among J. sabina populations across the vast region from Europe to Kyrgyzstan appear to be primarily shaped by geographical factors rather than taxonomic distinctions.Consequently, this leads to a relatively weak geographical pattern in the observed morphological variation [22].

Population History of J. sabina Varieties
The cpDNA-based analyses revealed that the J. sabina range was populated by different phylogenetic lineages corresponding to different varieties: var.sabina in the Cantabrians (SP7 population), Alps (AU3, AU4, IT1), Carpathians (RO3, PL), Crimean Mts.(CR) and Caucasus (GEO), var.balkanensis in the Apennines (IT3, IT4), Carpathians (RO1), Apuseni Mts.(RO2), and Balkans (CRO, NM, GR, BG), as well as the Asian variety in Tian-Shan (KYR) ([14] this study).These lineages formed three haplogroups in the MSN in the present study: one haplogroup involved 12 haplotypes of var.sabina from Europe and the GEO population, the second comprised two Kyrgyz haplotypes, and the third included five haplotypes of tetraploid var.balkanensis.As the haplogroups differed from one another by numerous mutational steps, the haplotype network was congruent with category I of Avise [55], i.e., the gene tree of allopatric lineages distinguished by prominent genetic gaps.This finding was supported by the multimodal mismatch distribution for the total dataset and bimodal distribution for the set containing diploid populations from Europe and Asia.Although multimodal/bimodal mismatch distributions can result from constant population sizes, they can also indicate the presence of different phylogenetic lineages [56].The European and GEO diploids underwent a population expansion, evidenced by unimodal mismatch distribution and significantly negative Tajima's D and Fu's Fs parameters.All these localities shared the H3 cpDNA haplotype, which indicates their common ancestry.The highest haplotype (Hd = 0.86) and nucleotide (π = 0.00048) diversities of the GEO population [14] can imply that it is older than the European diploid localities.GEO could not exist in the middle/late Miocene when the present-day areas of Georgia and Azerbaijan were covered by the Parathetys Sea; however, the climate cooling accompanied by the sea regressions promoted xerophytic vegetation in Eurasia in the late Miocene [57,58].
Inference of the time and place of var.balkanensis origin is difficult because of its hybrid ancestry followed by polyploidisation [11].As the ranges of var.balkanensis and its both parental taxa, i.e., maternal diploid var.sabina and paternal tetraploid J. thurifera-like species, do not overlap currently, it was hypothesised that the interspecific hybridisation was ancient [11].The fossils of potential J. thurifera paternal progenitor, Juniperus pauli Z. Kvaček, were found in the present Czechia and dated to the transition of Eocene/Oligocene [59].Juniperus pauli grew in humid, warm-temperate forests, while modern J. thurifera populates mountain regions with low precipitations and substantial sun radiation [14].The size stability of the var.balkanensis genome implied a long-lasting 'triploid bridge' between ancient hybridisation and recent tetraploidisation, as remarkable downsizing is expected after ancient polyploidisation [8].The SNP-and SilicoDArT-based hierarchical clustering revealed that the GR population split off from a common ancestor first among the studied var.balkanensis populations [14].Except for GR, populations of var.balkanensis in the Balkan and Apennine Peninsulas shared haplotype H14, the most common and centrally located in the tetraploid haplogroup.Thus, this haplotype could have originated in a single location and spread throughout Southern Europe.Population expansion in the Balkan and Apennine Peninsulas was evidenced by unimodal mismatch distribution of cpDNA sequences and significantly negative Tajima's D and Fu's Fs parameters.The exchange of flora between the eastern, western and southern parts of the Mediterranean Basin was facilitated at the end of the late Miocene due to the almost complete drying up of the sea, referred to as the Messinian Salinity Crisis [60].The refilling of the Mediterranean basin 5.33 Mya reduced plant migrations [61].

Conclusions
Within the area from Spain to Kyrgyzstan, we identified three allopatric phylogenetic lineages of J. sabina using cpDNA haplotypes.These lineages are congruent with three varieties: diploid var.sabina found in Europe and Georgia, tetraploid var.balkanensis from the Apennine and Balkan Peninsulas, and a distinct diploid variety in Kyrgyzstan [11][12][13][14].European var.sabina and var.balkanensis underwent population expansion.Although the juniper varieties are distinct based on the cpDNA markers, the population structure inferred from the nuclear microsatellites was unclear, since some genetic clusters included geographically distant populations or populations with different ploidy levels.Nevertheless, this outcome was congruent with the previous population structure assessment based on the SNP and SilicoDArT markers [14] and morphological parameters [22].We suppose that clonal propagation, inbreeding and limited gene flow between populations occupying different mountain ranges contribute to a moderate genetic variation of the J. sabina populations.The most intriguing result of our study was the finding that var.balkanensis occurred on the K-rich soils.The ample K resources in the soil probably alleviate the detrimental impacts of drought, thus enabling var.balkanensis to develop under hot and dry climatic conditions.

Figure 3 .
Figure 3. Genetic structure of 16 J.sabina populations (K = 5) inferred by Puechmaille's approach[36] with nuclear microsatellite markers' dataset.Each colour corresponds to one of five genetic groups.

Figure 4 .
Figure 4.The minimum spanning network (MSN) of 19 cpDNA (H1-H19) haplotypes in the J. sabina populations.Haplotypes were established by[14].CR-Crimea; the remaining population names are according to Table1.Blue colour indicates haplotypes of var.sabina, red-var.balkanensis, and green-the Asian variety from Kyrgyzstan.

Figure 4 .
Figure 4.The minimum spanning network (MSN) of 19 cpDNA (H1-H19) haplotypes in the J. sabina populations.Haplotypes were established by[14].CR-Crimea; the remaining population names are according to Table1.Blue colour indicates haplotypes of var.sabina, red-var.balkanensis, and green-the Asian variety from Kyrgyzstan.

Figure 5 .
Figure 5. Mismatch distributions of the cpDNA sequences of J. sabina based on pairwise nucle differences in the total sample (A) and three subsamples: all diploids (B), European diploid GEO (C), and tetraploids (D).Solid lines indicate the observed distributions under the popul growth-decline model and dotted lines indicate the expected distributions under this model.

Figure 5 .
Figure 5. Mismatch distributions of the cpDNA sequences of J. sabina based on pairwise nucleotide differences in the total sample (A) and three subsamples: all diploids (B), European diploids and GEO (C), and tetraploids (D).Solid lines indicate the observed distributions under the population growth-decline model and dotted lines indicate the expected distributions under this model.

Figure 6 .
Figure 6.Box plots showing differences in soil parameters between areas inhabited by the sa and balkanensis (B) varieties of J. sabina in Europe: N-nitrogen, P-phosphorus, K-pot CaCO3, EC-electrical conductivity, OC-organic carbon, pH; md-median differen

Figure 6 .
Figure 6.Box plots showing differences in soil parameters between areas inhabited by the sabina (S) and balkanensis (B) varieties of J. sabina in Europe: N-nitrogen, P-phosphorus, K-potassium, CaCO 3 , EC-electrical conductivity, OC-organic carbon, pH; md-median difference.The horizontal lines in the boxes show the median values.The ends of the boxes are on the first and third quartiles.The upper and lower whiskers show the maximum and minimum values, respectively.
horizontal lines in the boxes show the median values.The ends of the boxes are on the first and third quartiles.The upper and lower whiskers show the maximum and minimum values, respectively.

1 .
Chemical Soil Differences in the European Range of J. sabina

HWE) Lat (N) Lon (W/E) EUROPE
# Leaf samples were taken from the Kórnik Arboretum of the Institute of Dendrology of the Polish Academy of Sciences.

Table 2 .
Diversity parameters of three nuclear microsatellite loci studied in J. sabina.N A -number of alleles, H O -observed heterozygosity, H S -gene diversity, G IS -inbreeding coefficient, P(HWE)-pvalue for HWE.

Table 3 .
Heat map for pairwise genetic differentiation (Rho; ρ) values between the J. sabina populations based on the nuclear microsatellite markers.*-values statistically significant after the sequential Bonferroni's correction.A deeper shade of red indicates a higher genetic differentiation.

Table 4 .
Demographic histories of the total sample as well as diploid and tetraploid populations of J. sabina inferred from cpDNA markers by different statistic parameters (sp) under the assumption of neutrality.inds-individuals.