Three systems of molecular markers reveal genetic differences between varieties sabina and balkanensis in the Juniperus sabina L. range

Abstract

Key message

J. sabina var. balkanensis, that is of hybrid origin, and its maternal progenitor J. sabina var. sabina are genetically distinct with respect to cpDNA, SNP, and SilicoDArT loci. Mostly non-overlapping distributions of the sabina and balkanensis varieties are the result of their different climatic requirements.

Context

Juniperus sabina L. is present in the Eurasian mountains, but its range is severely fragmented. In Europe, two varieties of J. sabina occur: var. sabina and var. balkanensis, the latter being an allotetraploid hybrid between the diploid var. sabina and a tetraploid ancestor of Juniperus thurifera L. The distribution of the two varieties is mostly disjunct.

Aims

Assess the taxonomic affiliation and genetic differentiation of the populations of var. sabina and var. balkanensis in Europe and Asia using cpDNA, SilicoDArT, and SNP markers. Identify climatic niches of both juniper varieties in Europe.

Methods

Altogether, 21,134 SilicoDArT, 8,579 SNP, and four cpDNA loci were used. Seven climatic variables were compared in sites inhabited by var. balkanensis and the two parental species.

Results

The SilicoDArTs and SNPs revealed a pattern of population differentiation that was congruent with the cpDNA analysis. The hybrid var. balkanensis occupies habitats with higher temperatures and intermediate levels of precipitation compared to both parental taxa.

Conclusion

The low genetic variation and significant genetic differentiation among J. sabina populations likely result from the restriction of gene flow imposed by the mountain ranges. The balkanensis variety is able to cope with hot and dry climates probably thanks to the admixture of J. thurifera genes.

1 Introduction

Climate changes, a recurring phenomenon throughout Earth’s history, are the primary factors driving species extirpation, diversification, redistribution, and the subsequent structuring of genetic variation within species' ranges (Arenas et al. 2012). The direct result of population extirpation and species range fragmentation is an increase in inter-population distances, which as a consequence impede or even interrupt regional gene flow. When gene exchange is limited, genetic drift can result in the loss of genetic variation, potentially reducing the population's adaptive ability (Dobeš et al. 2017). On the other hand, new mutations and allele frequency changes resulting from the selection and genetic drift drive differentiation of the local populations. As a consequence of these processes, an isolation-by-distance pattern (IBD; Wright 1943), that is revealed by the positive correlation between geographic and genetic distances matrices, can be produced in heterogenous landscapes. In turn, the redistribution of species ranges enforced by climate change can result in the formation of hybrid populations, where individuals from different genetic lineages interbreed and produce mixed progeny (Barton and Hewitt 1989).

The species population histories are often very complex; thus, it is essential to utilize different molecular marker systems to provide reliable data on the extent and nature of genetic differentiation between the studied populations. Application of biparentally inherited nuclear markers and uniparentally inherited chloroplast (cp) or mitochondrial (mt) DNA fragments, which differ in respect to both mutation rate and effective population size, can explain the role of gene flow and selection in the species diversification, the emergence of adaptations, or the maintenance of species distinctness in the face of hybridization and introgression (Rieseberg et al. 1996). The modern genome-wide sequencing and bioinformatics methods allow to study the complexity of genomes and the role of evolutionary processes in shaping the genetic diversity of populations and species on an unprecedented scale so far (Adhikari et al. 2017).

Among terrestrial ecosystems, deserts and xeric shrublands are the most threatened by the present-day climate change, as warming and precipitation decrease can push them to their ecological limits (Tielbörger and Salguero-Gómez 2014; Li et al. 2018).

One of the shrub species adapted to arid and semi-arid biomes is Juniperus sabina L. (savin juniper), being widely distributed in Eurasia, but its range is strongly fragmented (Farjon 2005; Adams 2014). This juniper is a dioecious, broadly prostrate, sometimes decumbent shrub, reaching up to 1 m in height, and able to cover an area of about 10–12 m in diameter or even more (Farjon 2005). This species exhibits intensive clonal dissemination, with individual shrubs often outcompeting each other. Populations are scattered and typically confined to several hectares in size. According to Adams (2014), two varieties can be distinguished in the J. sabina range: sabina and balkanensis R.P. Adams and A.N. Tashev, although they hardly differ morphologically and are difficult to recognize in the field (Mazur et al. 2021a). The var. balkanensis is allotetraploid, likely as a result of the ancient hybridization between var. sabina and ancestor of Juniperus thurifera L., which are the maternal diploid and paternal tetraploid parents, respectively (Adams et al. 2016). The presence of var. balkanensis has been reported in Italy, the Balkan Peninsula, and western Turkey (Adams et al. 2016, 2017, 2018; Rajčević et al. 2022). Individuals of var. balkanensis were also found in one population in Spain, where they coexist with var. sabina (Farhat et al. 2020a). Today, the two parental species of var. balkanensis are sympatric only in the mountains of the Iberian Peninsula and in the Western Alps, while further east only var. sabina occurs (Farjon 2005; Adams 2014). J. thurifera probably occurred in a wider range in the past (Adams et al. 2016), but today the species can be found in mountains and within loose forests in the western Mediterranean region. The European range of J. thurifera includes isolated patches in Spain and southern France. The species is heliophilous and xerothermic, growing under continental and cold Mediterranean climates (Farjon 2005; Adams 2014). J. sabina is a mountainous species, inhabiting coniferous forests and invading meadows and pastures, and it can also be found on rocky mountain slopes. The species prefers abundant sunlight and is adapted to a climate with dry and warm summers and relatively cold winters (Farjon 2005; Adams 2014).

In this study, three systems of molecular markers were used to investigate the genetic diversity and differentiation among populations covering the majority of range of J. sabina. These were: cpDNA, SNP, and SilicoDArT markers. SNPs and SilicoDArTs were obtained through the DArTseq technology, combining DArT (Diversity Array Technology) with NGS (next-generation sequencing). The DArT method was developed by Diversity Array Technology Pty Ltd. (DArT, Canberra, ACT, Australia) for whole-genome profiling without any prior sequence information (Jaccoud et al. 2001; Kilian et al. 2012). The present study aimed to (1) identify taxonomic affiliation of the J. sabina populations from the Romanian mountains and Asian populations from Kyrgyzstan and Georgia based on the cpDNA haplotypes, (2) assess if var. balkanensis and var. sabina differ in respect of the DArTseq markers, and (3) to identify the climatic niches of var. sabina and var. balkanensis in Europe.

2 Materials and methods

2.1 Sampling and DNA isolation

In total, 108 samples of J. sabina representing 17 populations were used for cpDNA analyses, and 94 specimens from 14 populations were used for the DArTseq investigation (Appendix Table 1; Fig. 1). The populations are located in different mountains: Cantabrians, Alps, Apennines, Balkans, Apuseni, Carpathians, Crimean, Caucasus, and Tian-Shan. Samples were collected from individuals clearly delimited in the field, usually growing no closer than 30 m from each other. The populations from Spain (SP7), Austria (AU3, AU4), Italy (IT1), Ukraine (CR), and Poland (PL) were previously identified as var. sabina, whereas Italian (IT3, IT4), Bulgarian (BG), Croatian (CRO), Greek (GR), and North Macedonian (NM) as var. balkanensis, based on the cpDNA markers (Adams et al. 2017, 2018). Populations from Romania (RO1, RO2, RO3) and those from Kyrgyzstan (KYR) and Georgia (GEO) have not been studied so far with molecular markers.

Fig. 1

Locations of J. sabina var. sabina (open circles) and J. sabina var. balkanensis (black circles) populations in Europe and Asia. Population codes according to Appendix Table 1

To obtain genomic DNA, the top of the leafy twig was taken from each studied tree. Twigs were preserved in plastic bags with silica gel and stored at room temperature until DNA extraction. Before DNA isolation, the plant tissues were homogenized with a TissueLyser mill (Qiagen, Hilden, Germany). DNA extraction was carried out using a Qiagen DNeasy Mini Kit (Qiagen) following the manufacturer’s instructions. The DNA concentration in each sample was established between 50 and 100 ng μL⁻¹ with the help of a NanoDrop 2000 spectrophotometer (Thermo Scientific, Waltham, MA, USA).

2.2 Molecular analyses

2.2.1 Chloroplast DNA analyses

For cpDNA analyses, four intergenic spacers, trnL-trnF, trnD-trnT, trnS-trnG, and petN-psbM (Appendix Table 2), were used. Proportions of the PCR components and amplification program were described by Mazur et al. (2021b). The cpDNA regions were amplified in a Labcycler Basic (SensoQuest, Göttingen, Germany). PCR products were cleaned up with an EPPiC Fast enzymatic mixture (A&A Biotechnology, Gdańsk, Poland) and then sequenced using a BrilliantDye Terminator v. 3.1 Cycle Sequencing kit (NimaGen, Nijmegen, The Netherlands). An ExTerminator Kit (A&A Biotechnology) was used to purify the DNA sequencing reactions. The order of the nucleotides in cpDNA fragments was determined with an ABI PRISM 3500 Genetic Analyser (Thermo Scientific, Waltham, MA, USA). Variable sequences for different loci were deposited in GenBank with the accession numbers listed in Table 3 in the Appendix.

2.2.2 Genome-wide genotyping and data filtering

DNA samples of 94 individuals were sent to Diversity Arrays Technology Pty Ltd. in Australia (https://www.diversityarrays.com) for the sequencing and identification of genome-wide markers using the DArTseq method. In this method, genome complexity reduction is conducted by applying a combination of restriction enzymes to find and separate the most polymorphic genome regions, which are usually low-copy sequences (Kilian et al. 2012). To detect reproducible DNA polymorphism in J. sabina, the combination of PstI and MseI was chosen. The ends of digested fragments created by PstI were ligated to adapters and amplified by PCR and then sequenced on NGS platforms (Kilian et al. 2012; Andrews et al. 2016).

All DNA samples of J. sabina sent for DArTseq analyses showed good quality, and we received two types of markers that are available at the Zenodo repository (Jadwiszczak et al. 2023). These were SilicoDArTs, a type of dominant marker, and single nucleotide polymorphisms (SNPs), which are codominant marker systems. In the SilicoDArT loci, markers are scored ‘1’ for presence or ‘0’ for absence, and the marker is assigned as ‘- ‘ in the case when counts are too low to score confidently as ‘1’. In the SNP loci, the reference allele is assigned as ‘1’, and the SNP allele as ‘0’. The initial SilicoDArT data file comprised 241,216 markers, and the SNP data file included 100,346 markers. Next, data filtering was conducted using the R package dartR version 4.0.2 (Gruber et al. 2018) to reject poor-quality markers. The quality of both marker types was verified based on their reproducibility of 100% and call rate at the threshold of 95% for the SilicoDArTs and 85% for SNPs. Additionally, one ratio parameter at the threshold of 0.05 was also taken into account in the case of the SilicoDArT dataset. Reproducibility specifies the proportion of technical replicate assay pairs for which the marker score is consistent. The call rate is the proportion of samples for which the genotype call is either ‘1’ or ‘0’. One ratio determines the proportion of samples for which the genotype score is ‘1’.

2.3 Climatic data

For four sites of var. sabina, four var. balkanensis, and four localities of J. thurifera (Appendix Table 4), seven climatic parameters were analyzed. The parameters were as follows: mean annual dew/frost point temperature (°C) at 2 m above the earth's surface (DFP), maximum (Tmax) and minimum (Tmin) air temperatures (°C) at 2 m above the earth’s surface, mean annual soil humidity within the root zone (the layer from the surface 0 cm to 100 cm below grade; 0%—completely water-free soil, 100%—completely saturated soil, RZW), mean annual cloudiness (%; CA), mean daily amount of precipitation at the soil surface over the year (mm day⁻¹; PC), and the mean daily photosynthetically active radiation at the ground surface over the year (W m⁻²; PAR). All climatic data came from the POWER Release 901 project of the National Aeronautics and Space Administration (NASA) Langley Research Center (https://power.larc.nasa.gov/data-access-viewer/; accessed 13 November 2022). The data are based upon the Modern Era Retrospective-Analysis for Research and Applications (MERRA-2) assimilation model (Bosilovich 2008). The MERRA model is a general circulation model that shapes estimates of atmospheric variables based on assimilated and optimized observational data (Westberg et al. 2013). NASA provides the POWER data on a global grid with a spatial resolution of 1.0° latitude by 1.0° longitude for the solar radiation parameter, and with a resolution of 0.5° latitude by 0.625° longitude for the six remaining climatic parameters considered in this study.

2.4 Data analyses

2.4.1 CpDNA diversity and phylogenetic analysis

CpDNA haplotypes of 108 J. sabina specimens, including both polymorphic sites (parsimony informative sites and singletons) and indels, were established with the help of DnaSP 6.12.03 (Rozas et al. 2017). The haplotype diversity (H_d) was calculated manually according to the formula of Nei (1987):

$${\mathrm H}_{\mathrm d}=\left(1-\sum{\mathrm f}_{\mathrm h}{^2}\right)\times\frac n{n-1}$$

where f_h is the frequency of haplotype h in the sample and n is the total number of individuals. Nucleotide diversity (π) was established in MEGA version 11 (Tamura et al. 2021) using the partial deletion option that allowed us to include in the calculations all sites with coverage ≥ 5%. Phylogenetic analysis of the cpDNA haplotypes was described in the Appendix.

2.4.2 Polymorphism of DArTseq markers, population differentiation, and structuring

Based on the filtered sets of both types of DArTseq markers, the frequency distributions of polymorphism information content (PIC) values were established. The observed heterozygosity (H_O) and gene diversity (H_S) were estimated for those J. sabina populations that included at least 5 individuals. These two indices, as well as the expected heterozygosity of the total population (H_T) and inbreeding coefficient (F_IS), were calculated for the SilicoDArT and SNP datasets. For both SilicoDArTs and SNPs, based on the Euclidean genetic distances, the genetic relationships between all juniper individuals were studied using principal coordinate analysis (PCoA) and the hierarchical clustering method utilizing Ward’s approach (Ward.D2; Ward 1963). Considering populations with at least 5 specimens, the genetic differentiation was measured with pairwise F_ST estimates, and sequential Bonferroni’s correction (Rice 1989) was applied to avoid spurious positives. F_ST values were also established between groups of populations occupying different mountain ranges (Alps, Carpathians, Apennines, Balkans, Tian-Shan) as well as between the European var. sabina and var. balkanensis and between them and the Kyrgyz population from Tian-Shan. All of the above statistical analyses were performed using the dartR package (Gruber et al. 2018). The number of migrants (Nm) among the mountain ranges and the European and Asian savin junipers was assessed using the formula of Slatkin and Voelm (1991):

$$\mathrm{Nm}=\frac{\left(\mathrm d-1\right)\left(1-{\mathrm{F}}_{\mathrm{ST}}\right)}{4\mathrm d{\mathrm{F}}_{\mathrm{ST}}}$$

where N means the effective size of each population, m is the proportion of migrants arriving in a population per generation, and d indicates the number of populations analyzed. To investigate the correlations between F_ST and geographic distance (km) matrices (isolation by distance, IBD), the Mantel test (Mantel 1967) was performed in PAleontological STatistics (PAST) version 4.11 software (Hammer et al. 2001).

For both SNP and SilicoDArT markers, the most likely number of independent genetic groups (K) was assessed using the STRUCTURE 2.3.4 software (Pritchard et al. 2000) with 10 independent runs for K values ranging from 1 to 14. After the initial testing of different run lengths, the length involving 100,000 iterations after 10,000 burn-in periods was chosen. The locality information (LOCPRIOR model) as well as the admixture model with correlated allele frequencies were utilised in the analyses. Additionally, the recessive model was implemented for the SilicoDArTs. The CLUMPAK software was used to sum up and show graphical representations of the STRUCTURE results (Kopelman et al. 2015). The most probable number of genetic clusters was chosen based on the ad hoc statistic ΔK (Evanno et al. 2005). Based on K = 2 for the SilicoDArTs, the mean value of the cluster membership coefficient (Q) was calculated for each population across 10 runs.

2.4.3 Climatic niche analysis

In each site, the mean annual values of five climatic parameters recorded within the period 1981–2021 were standardized according to the Z score method. The maximum and minimum values of temperature in each site and each year were also considered and standardized. In the PAST software, the standardized parameter values were used to conduct redundancy analysis (RDA) with 9,999 permutations to study the potential impacts of the climatic parameters (explanatory variables) on the occurrence of var. sabina, var. balkanensis, and J. thurifera (response variables).

3 Results

3.1 Genetic diversity, population differentiation, and structuring

3.1.1 CpDNA markers

The nucleotide sequence lengths were 675 bp, 649 bp, 822 bp, and 878 bp for the trnL-trnF, trnD-trnT, trnS-trnG, and petN-psbM fragments, respectively. The total length of the combined cpDNA sequence was 3,024 bp. Altogether, 33 nucleotide substitutions (8 singletons and 25 parsimony informative) and 18 indel positions were found. The length of indels ranged from 1 to 20 bp. Out of 19, 12 haplotypes (H1–H12) were revealed in diploid European and GEO localities, five (H13–H17) in the tetraploid populations, and two (H18, H19) in the Kyrgyz population (Appendix Table 5). The highest number of haplotypes (5) was found in the GEO population, which resulted in the highest H_d (0.86) and π (0.00048). The total sample was characterized by H_d = 0.82 and π = 0.00404.

The BI phylogenetic tree revealed two main clusters of haplotypes. One cluster with 100% bootstrap support included five haplotypes revealed in var. balkanensis populations that grouped with those of J. thurifera and J. sabina × J. thurifera triploid hybrids (Appendix Fig. 3). The second cluster with 100% bootstrap support was divided into the var. sabina haplotypes from Europe, Georgia, and Azerbaijan, and another group comprising two Kyrgyz haplotypes. The Iranian sample seemed to be an outgroup with respect to all accessions in the cpDNA tree.

3.1.2 DArTseq markers

As a result of filtering, 21,134 SilicoDArT and 8579 SNP markers were retained for the J. sabina specimens. Most of the SilicoDArTs (99.99%) showed PIC values that were above 0.1 (Appendix Fig. 4), with a mean value of 0.195. In turn, most of the SNPs (76.78%) showed PIC values below 0.1, with a mean PIC of 0.074.

At the population level as well as for the total sample, the values of H_O and H_S were higher for SilicoDArT than for SNP markers (Appendix Table 6). The highest H_O (0.1620) and H_S (0.1382) values for the SilicoDArTs were revealed in the North Macedonian population. The lowest values of these parameters were found in Austrian AU3, 0.0462 and 0.0372, respectively. For the SNPs, H_O values ranged from 0.0110 in AU3 to 0.0579 in RO1 from Romania, while H_S ranged from 0.0196 in AU3 to 0.0843 in Croatia. In general, H_O and H_S values were clearly higher for the tetraploid localities when compared to the diploid populations for both types of DArTseq markers. In the total sample, the average values of H_O and H_S were 0.1160 and 0.0977, respectively, for SilicoDArTs and 0.0393 and 0.0607, respectively, for SNPs. The value of the H_T parameter was also higher for SilicoDArT (0.1052) than for SNP (0.0749) markers. The inbreeding coefficient was high and positive (F_IS = 0.3611) for SNP loci and negative (F_IS =  − 0.18) for SilicoDArTs.

In the SilicoDArTs-based PCoA, a clear division of the two main groups was revealed. The first and more diverse group consisted of the individuals deriving from var. balkanensis populations, and the diploid specimens were in the second group (Appendix Fig. 5a). The group consisting of the diploid specimens included the clearly distinct Kyrgyz population from Asia. In this PCoA, the first axis explained 6.25%, and the second explained 2.79% of the overall variance. In PCoA based on the SNP markers, three distinct groups were revealed (Appendix Fig. 5b). The first group included individuals from the tetraploid populations, the second group comprised specimens from diploid European populations, and the third group was constituted by the Kyrgyz specimens. The first and second axes showed 9.46% and 5.24% of the total variance, respectively, in the PCoA distribution based on the SNPs.

In the hierarchical dendrogram based on the SilicoDArTs, the first main group included all diploid specimens, and the second group consisted of tetraploids. In the cluster with diploids, European specimens were separated from the Kyrgyz samples (Appendix Fig. 6a). In the SNP-based hierarchical clustering, two main groups were revealed: one with the Kyrgyz specimens and the second with the European specimens, which were divided into diploids and tetraploids (Appendix Fig. 6b).

For the SNP markers, all pairwise population comparisons were statistically significant (F_ST from 0.046 to 0.512; Appendix Table 7). In the case of SilicoDArT loci, F_ST values ranged from 0 to 0.188. Based on the SNPs, populations from different mountain ranges were significantly genetically differentiated, and the number of migrants was the lowest between Tian-Shan and all remaining groups (Nm = 0.1–0.3; Appendix Table 8). For SilicoDArTs, the highest F_ST values were noted between Tian-Shan and var. sabina populations from the Alps and Carpathians. Considering three groups of populations, i.e., var. sabina in Europe, var. balkanensis and the Kyrgyz population, both kinds of DArTseq markers revealed higher genetic differentiation and a lower number of migrants between Kyrgyzstan and var. sabina compared to Kyrgyzstan-var. balkanensis (Appendix Table 8). Isolation by distance (IBD) was detected for both SNP (r = 0.93, P = 0.004) and SilicoDArT (r = 0.83, P = 0.010) markers.

The most likely number of genetic clusters, based on the method designed by Evanno et al. (2005) was K = 4 and K = 5 for the SilicoDArTs and SNPs, respectively (Appendix Fig. 7). In the case of K = 4 for SilicoDArTs, the first group involved var. sabina individuals from European populations as well as the Kyrgyz specimens (green color), and the second group consisted of the var. balkanensis populations from Romania, Italy, Croatia, and North Macedonia (orange), and third and fourth groups were represented by Greek and Bulgarian specimens (different proportions of violet and orange colors), respectively. When K = 2 was considered for SilicoDArTs (Appendix Fig. 7a), the populations were divided into var. sabina and var. balkanensis. The mean values of cluster 2 membership coefficient indicated that var. balkanensis populations likely exhibited an admixture of J. thurifera genes, ranging from 16% (Italy) to 20% (Greece) (Appendix Table 9). Based on the SNPs with K = 5, the STRUCTURE assigned most Romanian and Italian individuals to one cluster (orange), some Greek and one specimen from North Macedonia to a second cluster (violet), while the remaining var. sabina and var. balkanensis individuals showed an admixture of at least two gene pools (Appendix Fig. 7b).

3.2 Climatic niches

The RDA clearly differentiated var. sabina, var. balkanensis, and J. thurifera sites (F = 13.21, p = 0.001; Fig. 2). The first axis in the RDA explained 49% of the variability in the relationship among the occurrence of the three junipers and climatic parameters. This axis separated J. thurifera sites from both var. sabina and var. balkanensis. The second axis was responsible for explaining 46.9% of the variability, and it differentiated var. sabina sites from var. balkanensis and J. thurifera sites. The occurrence of var. sabina revealed strong positive correlation with the precipitation (PC), cloudiness (CA), and soil humidity (RZW) parameters and negative correlations with the photosynthetically active radiation (PAR), minimum (Tmin) and maximum (Tmax) temperatures, and dew/frost point (DFP) variables. The very acute angles between PC, CA, and RZW indicated that they were strongly correlated with each other. The J. thurifera sites showed positive correlations with PAR, Tmin, and Tmax and negative correlations with the humidity parameters. The var. balkanensis sites were positively correlated with PC, CA, RZW, DFP, and Tmax.

Fig. 2

Redundancy analysis (RDA) plot showing the relationship between var. sabina (JSS1-4; squares), var. balkanensis (JSB1-4; dots), and J. thurifera (JT1-4; triangles) sites and the climatic variables in the studied sites. DFP dew/frost point temperature, Tmax and Tmin maximum and minimum air temperatures, respectively, RZW soil humidity within the root zone, CA cloudiness, PC precipitation, PAR photosynthetically active radiation

4 Discussion

4.1 Distribution and genetic differentiation of the J. sabina varieties

The previous investigations showed that var. sabina and var. balkanensis were clearly distinct with respect to the cpDNA markers (Adams et al. 2016, 2018; Farhat et al. 2020a). In the present study, these markers were used to identify J. sabina populations that have remained genetically indefinite thus far (Mazur et al. 2021a). In the Romanian RO1 from the Southern Carpathians and RO2 from the Apuseni Mts., the cpDNA haplotype H13 was found. This haplotype formed a cluster with other haplotypes found in var. balkanensis populations. Thus, the RO1 and RO2 populations belong to var. balkanensis, and RO2 is its northernmost population at present. The RO3 population located in the Eastern Carpathians shared the H3 and H4 haplotypes with other var. sabina populations from Europe. In the Georgian population, the H3 haplotype was revealed, being most frequent in var. sabina populations. The clustering of two haplotypes from Azerbaijan (Hojjati et al. 2018) with other var. sabina samples in the BI tree implied that the variety extends at least to the western side of the Caspian Sea. The H18 and H19 haplotypes from Kyrgyzstan formed a distinct subgroup within the diploid cluster of J. sabina, suggesting the existence of a distinct variety in this country. It is not excluded that the same Kyrgyz haplotypes were also found in Kazakhstan and Mongolia, as Adams et al. (2016) noted that the Kazakhstan/Mongolia cpDNA cluster was clearly distinct from that from Europe and Azerbaijan.

Examination of cpDNA haplotypes did not reveal a population where the two savin juniper varieties co-occurred. Instead, var. balkanensis was found in the Apennine and Balkan Peninsulas, while var. sabina was spread from Spain to Georgia. The sabina and balkanensis varieties coexist in Sierra de Baza in Spain (Farhat et al. 2020a). Haplotypes of var. balkanensis were nested with the GenBank accessions of J. thurifera and J. sabina × J. thurifera triploid hybrids (Farhat et al. 2020b) in the cpDNA tree. This is the consequence of the origin of var. balkanensis, which is the allotetraploid hybrid between diploid var. sabina and the tetraploid ancestor of J. thurifera (Adams et al. 2016; Farhat et al. 2019). The high similarity of the cpDNA fragments of J. thurifera and var. balkanensis indicates that J. thurifera must be the paternal parent, as chloroplasts are most likely transmitted via pollen in Cupressaceae species (Farhat et al. 2019).

In the hierarchical clustering and PCoA based on the SilicoDArT and SNP markers, var. balkanensis and var. sabina specimens from Europe were also clearly separated. However, the position of the Kyrgyz individuals depended on the marker used. Using SilicoDArTs, the Kyrgyz genotypes formed a subcluster in the diploid group. In the SNP-based analyses, the Kyrgyz individuals grouped independently to both types of European populations. In the STRUCTURE clustering based on the SNPs (K = 5) and SilicoDArTs (K = 4), the Kyrgyz population was classified with var. sabina from Europe. In general, each clustering for SilicoDArTs revealed very high contribution of var. sabina genes in all studied populations, but when K = 2 was considered var. balkanensis populations were characterized by up to 20% of J. thurifera admixture. This result was congruent with the hybridization scenarios between diploid var. sabina and tetraploid J. thurifera that assumed backcrossing of tetraploid hybrid progeny with diploid maternal plants (Farhat et al. 2019). As an effect, var. balkanensis clusters with var. sabina in the nuclear internal transcribed spacer (ITS) region analyses (Adams et al. 2016). For K = 4, the SilicoDArTs showed some level of differentiation within var. balkanensis group because the Greek and Bulgarian populations had the admixture of a unique gene pool. The cpDNA analysis showed that the Greek and Bulgarian populations shared the H16 haplotype that was absent in other var. balkanensis populations. We were not able to find a biological sense of any clustering (K = 2, 3, or 5) based on the SNPs because the suggested groupings did not reveal taxonomic or geographical patterns. These discrepancies between the SilicoDArT and SNP marker analyses likely resulted from the lower informativeness of SNPs (PIC = 0.074) compared to SilicoDArTs (PIC = 0.195). In DArTseq technology, the SNP locus is considered biallelic, with the most common nucleotide acting as the reference allele and the second nucleotide representing an alternate allele. Any failure in calling SNP loci is referred to as missing data (Gruber et al. 2018). Calling failure can appear when the SNP state is ambiguous by consensus (Gruber et al. 2018), which seems to be more likely in polyploid than diploid genomes. It was found that the SilicoDArTs from the diploid genomes were more useful than SNPs in the determination of putative progenitors of polyploid Aegilops L. and Triticum L. species (Edet et al. 2018).

The distribution of J. sabina is confined to the mountain ranges; thus, a restricted gene flow was expected among the populations over the species range. As it was presumed, the highest F_ST and lowest number of migrants (Nm = 0.1–0.3) were noted between the Kyrgyz population from Tian-Shan and all other mountain ranges. Based on both types of DArT markers, the lowest F_ST and highest Nm were revealed between var. balkanensis populations located in the Apennines, Carpathians, and Balkans. However, we suppose that this result can indicate shared ancestry rather than actual gene exchange among them. In general, F_ST values were higher for SNPs compared to SilicoDArTs but most of the F_ST values were statistically significant at both mountain group and population levels. Consequently, the significant isolation-by-distance (IBD) pattern was revealed by both kinds of markers. Significant IBD strongly implies spatially restricted gene flow due to the limitation of pollen and seed movements. This limitation can result from the presence of mountain ranges as geographical barriers. Meta-analysis of 179 plant species from Central Europe and the European Alps revealed significantly higher genetic differentiation among the populations inhabiting the mountains compared to those from the lowlands, and among the populations from the alpine zone compared to those from the subalpine areas (Reisch and Rosbakh 2021).

4.2 Climatic niches of the J. sabina varieties

The influence of climate on juniper distribution is evident, as their diversification during the Tertiary was linked to progressive climate drying (Farjon 2005). The current range of J. sabina spreads from Spain to Mongolia but it is strongly fragmented, and the varieties sabina and balkanensis show nearly completely nonoverlapping distributions. This geographical disjunction can result from different climatic requirements of var. sabina and var. balkanensis. Our RDA analysis revealed that var. sabina was present in areas with the greatest cloudiness and precipitation but the lowest temperatures. In contrast, var. balkanensis occurred in the sites with the highest maximum temperature and dew/frost point. With respect to the precipitation, var. balkanensis occupies the areas with lower humidity than that in the var. sabina sites. Analyses of bioclimatic factors in the Balkan Peninsula showed a substantial drought tolerance of var. balkanensis as two populations from the Balkan Mts. experienced at least 20% less precipitation than the sites located in the Dinaric Mts. (Rajčević et al. 2022). This tolerance seems to be an effect of J. thurifera genes because J. thurifera, that is the paternal ancestor of var. balkanensis (Adams et al. 2016), inhabits sites with the highest photosynthetic radiation and lowest precipitation.

The almost allopatric distribution of the two varieties likely dates back to at least the last Pleistocene glaciation. The balkanensis variety could have survived the glacial maximum in the refugia located in the Iberian, Apennine, and/or Balkan Peninsulas. However, during the Holocene warming, its northward spread was likely impeded by the harsh climate of the Pyrenees and the Alps. The sabina variety as a mountainous species can easily cope with relatively cold winters with abundant snowfall; thus, its populations could have survived the Pleistocene in different mountain ranges, including the periphery of the Alps (Schönswetter et al. 2005). After the glacial retreat, refugial populations of var. sabina recolonized the upper mountain elevations. In a consequence of the long-lasting geographical separation of var. sabina and var. balkanensis, being a result of their different climatic adaptations, the clear genetic differentiation in the cpDNA, SNP, and SilicoDArT loci was shown (Adams et al. 2016; this study). The long-term isolation caused the genetic differentiation of several mountainous plants. It was shown that the populations of Campanula alpina Jacq. inhabiting the Carpathians and the Alps as well as C. alpina populations from the different Carpathians ranges were genetically differentiated (Ronikier et al. 2008). Genetic differentiation was also found between the populations of Pinus mugo Turra from the Sudetes, Carpathians, and Alps (Boratyńska et al. 2014).

In the face of current climate deterioration, the conclusions regarding the further fate of J. sabina are uncertain. The species has survived many climatic changes thus far. On the one hand, based on SNPs and SilicoDArTs, var. sabina appears to be less genetically diverse, making it potentially more vulnerable to climate warming than the allotetraploid var. balkanensis. It was found that the hybrid populations having mixed gene combinations could be more resilient to climate change than the parental taxa (Brauer et al. 2023; Turbek and Taylor 2023). On the other hand, it was reported that var. sabina intensively overgrew pastures in Spain (García-Cervigón et al. 2017) and in the Alps, and it was even proposed to cut down bushes to stop this process (Reinhard 2017). In the Carpathians and the Dynaric Mts., the shepherds burn bushes of var. balkanensis to prevent the overgrowth of grassy pastures by the junipers (Łuczaj and Ronikier, pers. inf.). Climate change, but also the removal of shrubs by humans, can force both varieties of J. sabina to move to higher altitudes, to the subnival zone in the mountains. In our opinion, the juniper populations in the low mountain ranges are at the greatest risk.

4.3 Genetic resources

At the population level, the values of the diversity parameters H_O and H_S for both types of DArT markers were up to fivefold higher in var. balkanensis populations compared to var. sabina populations. The clear distinction of heterozygosity measures between the diploid and tetraploid junipers was rather expected, as previous molecular analyses revealed that 75% of var. balkanensis individuals and 47% of var. sabina were polymorphic with respect to the nuclear ITS regions (Adams et al. 2018; Farhat et al. 2020a). One of the most surprising results in our study was the quite low H_T values assessed using both SNP (0.0749) and SilicoDArT (0.1052) markers, although J. sabina individuals originated from 14 populations spread over a distance of 6600 km. The H_T parameter was similar or even higher in small-scale studies of other species: the palm Mademia argun (Mart.) Wurttenb. ex H.Wendl. from Sudan (0.036 and 0.127 for SNPs and SilicoDArTs, respectively; Elshibli and Korpelainen 2021) or the flowering tree Trema orientalis (L.) Blume in Uganda (0.40 and 0.06 for SNPs and SilicoDArTs, respectively; Nantongo et al. 2022). There are three possible explanations for this phenomenon. First, the low H_T value for SNP markers in J. sabina resulted from inbreeding, as F_IS was equal to 0.3611. Based on both SilicoDArT and SNP markers, the T. orientalis population was not inbred (Nantongo et al. 2022). Second, the low diversity level of the savin juniper could be a consequence of the lower PIC values of both markers compared to other investigations. The average PIC values were 0.17 for SNPs and 0.22 for SilicoDArTs in T. orientalis (Nantongo et al. 2022). Low PIC value of the DArT markers can result, among others, from the presence of null alleles, i.e., alleles that are not sequenced due to the mutation (Andrews et al. 2016). The third possible explanation is the much higher rate of molecular evolution per unit time in angiosperm than in gymnosperm species (De La Torre et al. 2017). Contrary to the variation of DArT markers, the cpDNA diversity was higher in var. sabina (H_d = 0–0.86, π = 0–0.00048) populations compared to var. balkanensis (H_d = 0–0.60, π = 0–0.00009). This phenomenon can likely be explained by the mixing of different phylogenetic lineages or survival in the isolated glacial refugia; however, this is a topic of our next study.

5 Conclusions

The SilicoDArTs and SNPs showed that the genetic variation in the savin juniper populations spread over a large area was low. Moreover, the parameters of genetic diversity were much lower in var. sabina populations compared to var. balkanensis. As genetic diversity is a prerequisite for populations to cope with environmental changes, the var. sabina populations, especially those from the lower mountain elevations, can be at a higher risk than var. balkanensis in the face of ongoing climate warming. The var. sabina inhabits sites with high precipitation and low temperatures; thus, it may suffer from rising temperatures and subsequent drying of arid and semiarid habitats in the future. As moisture excess during winter-spring season at high altitudes can limit the upward expansion of J. sabina (García-Cervigón et al. 2018), the species’ genetic resources should be preserved by ex-situ conservation measures. Successful ex-situ germination efforts were undertaken in the Kórnik Arboretum of the Institute of Dendrology of the Polish Academy of Sciences where var. sabina from the Pieniny Mts. is preserved.

Appendix

Table 1 Population names, geographical locations, and number of J. sabina individuals studied in respect of the DArTseq markers and cpDNA

Table 2 Sequences of chloroplast (cpDNA) primer pairs used in the J. sabina study (see Adams et al. 2016 and references therein)

Table 3 GenBank accessions of J. sabina, J. thurifera, and J. sabina × thurifera triploid hybrids used in the construction of phylogenetic tree from the Bayesian inference analysis using the concatenated cpDNA sequences

Table 4 Locations of J. sabina and J. thurifera populations analyzed in respect of seven climatic parameters

Table 5 Genetic variation of the concatenated trnL-trnF, trnD-trnT, trnS-trnG, and petN-psbM cpDNA fragments in 17 populations of J. sabina

Table 6 Diversity parameters of the SilicoDArT and SNP markers in the J. sabina populations

Table 7 Heat map for pairwise genetic differentiation (F_ST) values between the J. sabina populations based on SNP (above diagonal) and SilicoDArT (below diagonal) markers. Population codes according to Table 1

Table 8 Heat map for pairwise genetic differentiation (F_ST) values between the mountain ranges and varieties of J. sabina based on SNP (above diagonal) and SilicoDArT (below diagonal) markers

Table 9 The values of membership coefficient Q in the particular populations of J. sabina based on the K = 2 clustering using the SilicoDArT markers. Population codes according to Table 1

1.1 Phylogenetic analysis of the cpDNA haplotypes

Phylogenetic analysis was conducted in the program MrBayes 3.2.7a (Ronquist et al. 2012), which performs Bayesian inference (BI) using a variant of Markov chain Monte Carlo to approximate the posterior probabilities of the tree. In this analysis, the combined cpDNA fragments of three J. sabina individuals from Azerbaijan and Iran [14316 (BAYLU), 14317 (BAYLU), 101434 (TARI); Hojjati et al. 2018], six J. thurifera (JT2, JT3, JT5, JT9, JT16, JT39), and three J. sabina × J. thurifera triploid hybrids (JTS8, JTS14, JTS19; Farhat et al. 2020b) were also included (Appendix Table 3). The HKY substitution model (Hasegawa et al. 1985), established in MEGA based on the lowest AIC value, was used in the construction of the BI tree. In MrBayes, two independent runs of 10 million generations each were carried out. Each run consisted of three heated chains and one cold chain. Trees were sampled every 1000 generations. After excluding 25% of the first runs as burn-in, a majority-rule consensus tree with posterior probabilities (PPs) from two runs was obtained.

Fig. 3

Bayesian tree based on the concatenated trnL-trnF, trnD-trnT, trnS-trnG, and petN-psbM cpDNA fragments of J. sabina, J. thurifera, and J. sabina × J. thurifera triploid hybrids

Fig. 4

Frequency distribution of polymorphic information content (PIC) values for SilicoDArT and SNP markers retained after quality analyses in the set of 94 J. sabina individuals

Fig. 5

Principal Coordinates Analysis (PCoA) revealing the Euclidean genetic distances between 94 J. sabina individuals based on the SilicoDArT (a), and SNP (b) markers. Population codes according to Appendix Table 1

Fig. 6

Results of the hierarchical cluster analyses of 94 J. sabina individuals, utilizing Ward’s approach for SilicoDArT (a), and SNP (b) markers. Population codes according to Appendix Table 1

Fig. 7

Bayesian clustering of 94 J. sabina individuals from 14 populations based on a SilicoDArT loci with K = 2, K = 3, and K = 4, and b SNP loci with K = 2, K = 3, and K = 5. Population codes according to Appendix Table 1