Juniperus excelsa s. str. in crimea – differentiation and history inferred from genetic and morphological markers

Seven Crimean populations of J. excelsa were explored during the summer 2017, and additional material from the small, isolated and unstudied population in Chatyr Dag was obtained in 2019. In all populations, material from 59 trees was collected (Tab. 1, Fig. 1). Altogether, 245 cones and the same number of seeds, as well as 297 small parts of shoots with leaves were examined for their morphology. The leaves of 48 individuals secured in silica-gel were used for genetic tests, because the quality of remaining samples was low (Tab. 1).

Figure 1

The collected material was dried without pressing before measurement. The choice of traits analyzed (Tab. 2) was based on previous examinations by different authors, performed for J. excelsa and other species of juniper (Marcysiak et al. 2007; Mazur et al. 2010, 2018; Douaihy et al. 2012). Seven traits were examined under the Delta Optical SZ-630 microscope: cone length (CL), seed length (SL), seed width (SW), thickness of the last ramification shoot with leaves (ST), and two measurements of the width of the cone at its widest point, perpendicular to each other, the average of which gave the diameter of the cone (CD). Leaves were counted (LN) under the microscope on the last 5 mm long section of the shoot. The Delta Optical DLTCamViewer software was used to obtain digital images and to measure elements. The results were exported to html files, which in turn were copied to the database in Excel. Numbers of cone scales (CSN) and seeds (SN) were counted manually. Eight ratios were calculated with the measured/counted characteristics (Tab. 2).

The symmetry and unimodality of the distribution frequency of measured characteristics were checked using Shapiro-Wilk’s test, and mutual correlations of characteristics were examined using Spearman’s rank correlation coefficient (Zar 1999). The basic descriptive statistics (arithmetic means, minima and maxima, standard deviations, and coefficients of variations – CV) were obtained for all data and for each population. The set of individual means was created, as a basis for further analyzes, and the normality of traits distribution was checked again. ANOVA helped to define differences between populations, after verifying the homoscedasticity of variances with the Brown-Forsythe test. Differences between groups and populations were tested by Student t test or Tukey’s test, when characters fitted parametric tests, and U Mann-Whitney test or Kruskal-Wallis test, when they failed parametric tests; significance (P) was reported (Zar 1999). The influence of ecological and geographic factors on the differences between populations was investigated by comparing coastal and mountain populations as well as eastern and southernmost populations (Tab. 1). Discriminant analysis conducted for six populations (Chatyr Dag, Fiolent, Laspy, Nikita, Sudak, Uzundza) helped to verify the grouping of the populations. The Kara Dag and Bogatoye Ushelye populations, comprising only few cones, were excluded from these analyses. Traits used in this analysis had normal distributions and were chosen based on their significance in differentiating populations, showed by ANOVA, and the lack of strong correlations (Sokal and Rohlf, 2003). The software package STATISTICA 13.3 (StatSoft) was used for calculations.

An AX Plant Kit (A&A Biotechnology) was used for DNA isolation, according to the manufacturer’s instruction. To screen possible genetic variation, four non-coding cpDNA fragments [trnL-trnF (Taberlet et al. 1991), trnD-trnT (Grivet et al. 2001), trnS-trnG (Adams and Kauffmann 2010), and trnH-psbA (Shaw et al. 2005)] and the nuclear internal transcribed spacer region ITS1-5,8S-ITS2 (ITS) were sequenced. Amplifications of the chloroplast regions were conducted in 9 μl reaction mixture containing 4.0 μl of Multiplex PCR Master Mix (Qiagen), 2.2 μl of RNase-free water (Qiagen), 0.8 μl of primer mixture (0.2 μM of each primer), and 2.0 μl of diluted DNA. The total volume of the mixture for amplification of the ITS fragments was 15 µl [6.75 μl of Multiplex PCR Master Mix, 3.95 μl of RNase-free water, 1.3 μl of primer mixture (0.2 μM of each primer), and 3.0 μl of diluted DNA]. The PCR program for cpDNA fragments was according to Shaw at al. (2005); in turn, the PCR profile for the ITS amplification was described by Zalewska-Gałosz et al. (2010). After PCRs, all products were purified with the Exo-BAP kit (EURx), and cpDNA products were sequenced directly using a BrilliantDye Terminator v. 3.1 Cycle Sequencing Kit (Nimagen). The ITS fragments were first run on 1.5% agarose gels in 1 × TBE buffer (4.5 h, 60 A) and stained with SimplySafe (EURx). Then, the ITS bands were excised and purified using a Gel-Out kit (A&A Biotechnology). Before visualization on an ABI 3500 sequencer (Applied Biosystems), an ExTerminator Kit (A&A Biotechnology) was used to remove unincorporated dNTPs and primers from the sequenced products. Variable sequences for different loci were deposited with the accessions numbers SAMN19591084 – SAMN19591104 within the BioProject PRJNA735631 in GenBank (Supplementary Tab. 1).

Sequences were assembled using BioEdit v. 7.0.4 (Hall 1999) and adjusted manually, based on the Gen-Bank accessions LC420837, LC420838, LC420839, and LC420850 of J. excelsa for the alignment of the trnD-trnT, trnL-trnF, trnS-trnG, and ITS regions, respectively. The trnH-psbA fragments were aligned based on the accession MH566939 of Juniperus chinensis.

The DnaSP program, version 6.0 (Rozas et al. 2017), was used to indicate parsimony-informative sites (P_i) in each variable cpDNA fragment and calculate the haplotype diversity (H_d) of cpDNA markers in each population. To check whether the observed frequency distributions of the mutations in the haplotypes resulted from range expansion, Tajima’s D (Tajima 1989) as well as Fu’s Fs test (Fu 1997) and Fu and Li’s D* and F* statistics (Fu and Li 1993) were calculated using DnaSP. Range expansion is likely when Tajima’s D < 0, Fs is significantly negative, while D* and F* parameters are not significant (Fu 1997). Additionally, mismatch distribution analysis (distribution of the observed number of differences between pairs of haplotypes) was carried out in DnaSP to test potential current expansion in the J. excelsa populations. Arlequin, version 3.5.2.2 (Excoffier and Lischer, 2010), was used to construct the minimum spanning tree for cpDNA haplotypes. In the construction of this network, only parsimony-informative sites were used.

Phylogenetic analyses were conducted using a maximum likelihood (ML) method implemented in MEGA X version 10.1 (Kumar et al. 2018) and Bayesian inference (BI) in MrBayes 3.2.6 (Ronquist et al. 2012), based on the two data sets. The first set contained concatenated trnD-trnT/trnL-trnF/trnS-trnG fragments derived from 10 J. excelsa individuals representing different cpDNA haplotypes in the present study as well as seven accessions from GenBank with herbarium numbers: 14742(BAYLU), 13721(BAYLU), 14155(BAYLU), 8785(BAYLU), 9433(BAYLU), 13720(BAYLU) and 8786(BAYLU) (Supplementary Tab. 2), another set included the ITS sequences. The trnD-trnT, trnL-trnF, and trnS-trnG fragments were chosen because they were also amplified and showed polymorphism in some other J. excelsa specimens inhabiting different parts of the species range (Hojjati et al. 2018) and one J. chinensis (Mao et al. 2010). For the ITS tree constructions, nine sequences from this study as well as 12 J. excelsa (LC420820, LC420840, LC420855, LC420800, LC420850, LC420830, LC420810, LC421168, LC421145, LC421218, LC421203, LC421035), two J. communis (LC420860, LC420950), and two J. chinensis (LMH566939, FJ980279) sequences from GenBank were used (Supplementary Tab. 1, 2). The best-fit models of nucleotide substitution for ML analyses were based on the Akaike information criterion (AIC) and tested with MEGA. Tamura 3-parameter (T92; Tamura 1992) and Hasegawa-Kishino-Yano (HKY; Hasegawa et al. 1985) models were indicated for cpDNA and ITS fragments, respectively. The reliabilities of ML trees were inferred using bootstrapping with 1,000 replicates. To construct BI trees, the GTR substitution model with gamma-distributed rate variation across sites and a proportion of invariable sites was implemented. Two independent runs of 10 million generations each were carried out, with each run consisting of four chains (three heated, one cold). Trees were sampled every 1,000 generations. A majority-rule consensus tree with posterior probabilities (PPs) from two runs was produced after excluding the first 25% of runs as burn-in.

In the entire data set, unimodal frequency distributions of values were found for CL, CD, SL, SW, ST, CL/CD, CSN/CD, and slightly skewed for SL/SW, CD/SW, LN/ST, and CD/ST. The characteristics CSN, SN, LN, CD/SN, CL/SL were not normally distributed. In the set of all individual means, the frequency distributions were unimodal for all traits, and their variances were homoscedastic. The variation of the studied characteristics was moderate. Of all parameters, CD/SN showed the highest variation (Tab. 2), with CV > 25% in every population. The trait CL/CD was the most stable one, with CV values ranging from 2.61% in Nikita to 7.75 in Laspy (Supplementary Tab. 3). The examined characteristics were moderately correlated. The strongest correlations (r > 0.8) were found for CL and CD, SN and CD/SN, and LN and LN/ST. Strong relations, with r > 0.6, occurred between CL and CD/ST, ST and CD/ST, SW and SL/SW, and SL/SW and CD/SW. Leaf and shoot characteristics (LN and ST) were not correlated with features of cones and seeds (Supplementary Tab. 4).

The results of ANOVA performed on the set of individuals means revealed that CL, CD, SW, LN, ST, SL/SW, CD/ST, and LN/ST features with p < 0.05 were responsible for differentiating populations; however, according to Tukey’s test results, some of them differed only between single populations. Uzundza, Sudak, and Chatyr Dag differed most from other populations (Supplementary Tab. 5). Based on the results of Student’s t test performed for individuals means, the coastal and mountain groups of populations differed significantly in terms of: LN, SL/SW, and LN/ST. This test, when performed for eastern and southernmost populations, showed significant differences in CSN and LN/ST. Non-parametric tests, applied to all data, showed a greater number of significant differences (Tab. 2).

The discriminant analysis based on the traits: CL, CD, SW, and ST, showed grouping of populations: Fiolent, Laspy, and Nikita in the center of the graph, while the remaining three populations took the position outside this group (Fig. 2a). On the plot of all cases, one rather loose cloud was visible, but the outside located populations were more compact, and ellipses of the confidence area (95%) of Chatyr Dag and Uzundza only slightly overlapped each other (Fig. 2b). The first discriminant variable was most strongly connected with ST and CD, and therefore, the separateness of Chatyr Dag depended on these features. The second variable was determined mainly by SW, and the differentiation of Sudak from Uzundza was based on these traits.

Figure 2

Four cpDNA sequences were obtained for each J. excelsa individual, and the total length of the combined trnL-trnF/trnH-psbA/trnD-trnT/trnS-trnG sequence was 2,643 nucleotides. The trnL-trnF and trnS-trnG markers showed no polymorphisms; thus, 10 cpDNA haplotypes were established based on the variation of trnD-trnT and trnH-psbA fragments. In these fragments, 2 and 10 variable sites were found, respectively. Two sites in the trnD-trnT and five sites in trnH-psbA regions were parsimony-informative (P_i). The most frequent one was the H1 haplotype (72.9%), which was found in all populations (Tab. 3). The H2 haplotype (8.3%) occurred in three populations, and the H3 haplotype (4.2%) was noticed in two localities. The remaining seven haplotypes were observed only once. Haplotype diversity ranged from H_d = 0 in both Bogatoye Ushelye and Chatyr Dag to 0.9 in Fiolent. In the minimum spanning tree, the central position was occupied by the H1 haplotype (Supplementary Fig. 1). The remaining haplotypes differed from H1 maximally by three mutational steps.

Among the tests conducted to detect potential range expansion, only Fu’s Fs value was significantly negative (Fs = −2.813, P = 0.036), and Tajima’s D (−1.596; 0.1 > P > 0.05), Fu and Li’s D* (1.424; 0.1 > P > 0.05), and F* (0.430; P > 0.1) were nonsignificant. Mismatch analysis displayed a unimodal distribution (Supplementary Fig. 2), which can imply a recent demographic expansion.

The ML and BI phylogenetic trees were almost identical to one another for both ITS cpDNA sequence sets; thus, only the ML trees are shown and discussed. Phylogenetic analyses were based on 2,113 nucleotide sites of the concatenated trnD-trnT, trnL-trnF, and trnS-trnG cpDNA fragments. In general, analysis of the cpDNA regions showed that the J. excelsa haplotypes described in the Crimean Peninsula formed one group with J. excelsa from Greece, Bulgaria, Turkey, and Lebanon (Supplementary Fig. 3). In this group, only samples representing H2, H5, and H6 haplotypes, found in the Fiolent, Nikita, and Laspy localities, were distinct to some extent and formed two subgroups within the main cluster with 61 and 62% bootstrap supports; however, one subgroup included a sample from Greece.

It was difficult to obtain nuclear ITS sequences of good quality, and only nine fragments representing J. excelsa individuals from Kara Dag, Sudak, Fiolent, Nikita, Laspy, and Chatyr Dag populations were used (Supplementary Tab. 1). The total length of the ITS alignment was 616 nucleotides. In the ITS tree, all J. excelsa samples from Crimea clustered with the majority of individuals of the same species from other geographical locations (66% bootstrap support; Fig. 3). Two accessions of J. chinensis (MH566939 and FJ980279) and two J. communis (LC420860 and LC420950) were significantly distinct (99% bootstrap supports) from the J. excelsa ingroup.

Figure 3

Phylogenetic analyses of the cpDNA fragments revealed that the J. excelsa samples from the Crimean Peninsula were closely related to the individuals from other parts of the species range: Bulgaria, Greece, Turkey, and Lebanon. In general, this result was supported by the ITS phylogenetic investigation as the majority of J. excelsa individuals fell into one group, although the bootstrap value was not high. This observation confirms the statement of Adams et al. (2016) that J. excelsa is a genetically uniform species. A significant genetic similarity of the J. excelsa populations occupying different areas can result from a shared ancestry and indicate that the species range was more continuous in the past, with only recent fragmentation (Yena et al. 2005; Cordova 2007). According to Adams et al. (2014, 2016), only J. excelsa samples from Lebanon were different to some extent, and they formed a strongly supported subgroup in the phylogenetic tree based on the nuclear ITS and cpDNA markers. The molecular differences of the Lebanese specimens were supposed to be an effect of gene inflow from J. polycarpos into J. excelsa, as these two species co-occur in this region (Adams et al. 2016). However, one [14155 (BAYLU)] of the Lebanese accessions used by Adams et al. (2014, 2016) clustered with other J. excelsa samples in the analyses of Hojjati et al. (2018) as well as in the present study. In turn, Hojjati et al. (2018), who investigated genetic relationships between the Juniperus populations in Iran using nuclear and cpDNA loci, found that some Iranian samples classified as J. excelsa were related to J. polycarpos var. turcomanica, others to J. seravschanica. The authors concluded that there was no J. excelsa s. str. in Iran, which has previously been suggested by Adams et al. (2014). In our ITS tree, four Iranian accessions (LC421168, LC421203, LC421035, and LC421145), were placed with other J. excelsa samples, while the LC421218 individual was clearly different. Most probably, this unexpected similarity of some ITS sequences from Iran to other accessions of J. excelsa in the present work can result from the relatively short ITS fragments studied and the small number of variable sites.

Multivariate analyses pointed out similarities between Fiolent, Nikita, and Laspy, that is, populations located in the southern part of the Crimean coast, at about 50 m above sea level. These three populations took the central positions on the scatterplots, and their variability was quite high as their confidence ellipse covered a large part of the whole variability. Additionally, the H2 cpDNA haplotype, being the second most common haplotype, was only found in the Fiolent, Nikita, and Laspy localities. It is likely that the substantial morphological variation and the presence of H2 haplotype in three southernmost populations result from their common ancestry, indicating that they can be remnants of a huge continuous population in the past.

Although the analysis showed little morphological difference between the groups of mountain and coastal populations, two of the mountain populations, Chatyr Dag and Uzundza, differed by the greatest number of features from all the others. The Chatyr Dag population, with the smallest and elongated cones and seeds, grows in the most severe climate conditions among all Crimean localities of the Greek juniper, and is isolated from others. Therefore, its distinctiveness may be caused by climatic conditions, as the effects of climate on the growth and diversity of junipers have already been described (De Soto et al. 2014; Pellizzari et al. 2014; Mazur et al. 2016). The other possibility is the loss of variability as a result of the decline of this population, what can be enhanced by the lack of cpDNA haplotype diversity in Chatyr Dag. The Uzundza population grows in mountains, but in the milder climate, and the distances between this population and its neighbors are smaller, when compared to the Chatyr Dag population. Despite this, morphological differences between Uzundza and the closest populations in Laspy and Nikita were found. Its haplotype diversity is quite substantial, which results from the presence of unique haplotypes, H7 and H8. Again, the morphological distinctiveness of this population can be explained either by the effect of the climate, that seems to be here most favorable to the species, or by genetic drift acting on this isolated locality. The comparison between the populations from the eastern and southern ends of the range of J. excelsa in the Crimea revealed few differences, but Sudak being the largest among the eastern populations differed in several features from all the others. The population was also genetically distinguished by having two unique haplotypes, H9 and H10.

The minimum spanning network constructed based on the parsimony-informative sites found in the trnH-psbA and trnD-trnT cpDNA regions in the J. excelsa individuals from the Crimean Peninsula represented a typical starburst topology, suggested for an abundant species that expanded recently from a small number of individuals (Avise 2000). In such networks, the central position is occupied by an ancestral haplotype, which is the most common and widespread haplotype in the studied area; the remaining haplotypes are assumed to be recent derivatives connected to the central one mostly by short branches (Avise 2000). In the Crimean part of the J. excelsa range, the H1 haplotype was found in each locality and its contribution was the highest; thus, it seems to be ancestral. Except for H2 and H3 haplotypes, all others occurred once in different populations. Previous minimum spanning networks including the J. excelsa accessions involved also additional Juniperus species and varieties from Europe and Asia. However, the J. excelsa haplotypes differed from one another by one to three mutations (Adams et al. 2014, 2016). On the other hand, the results of analyses performed to test signatures of demographic expansion of J. excelsa in the Crimean Peninsula were contradictory. The distribution of the observed number of differences between pairs of haplotypes was unimodal, implying recent expansion. Significantly negative value of Fu’s Fs and non-significant values of Fu and Li’s D* and F* also seemed to confirm a recent expansion scenario. In turn, Tajima’s D value was non-significant, albeit negative, and this did not agree with the expansion model. We think that lack of concordance between different tests conducted to detect the demographic history of Crimean populations of J. excelsa can result from the low level of polymorphism in the cpDNA fragments studied. Similar outcomes and conclusions have been described for the two vitex species Vitex rotundifolia and V. trifolia in East Asia, both representing distributions under an exponential growth rate model, but Tajima’s D and Fu’s Fs tests results were non-significant (Sun et al. 2019).

The homogeneous character of J. excelsa s. str. within the whole range, confirmed previously by several authors and in our study, indicates the same origin of its many scattered populations. The traces of the Mediterranean flora in Crimea could be dated back to the Miocene-Pliocene (Didukh 1992; Cordova 2007), and it is assumed that the Crimean Peninsula could have served as a refuge for thermophilous species during the cold stages of the Pleistocene (Didukh 1992; Cameron et al. 2013). During the Pleistocene, in Crimea, the ranges of species were changing not only due to temperature fluctuations, but also because of changes in the level of the Black Sea, which remained low until the Holocene, what enabled plants migrations (Cordova 2007; Krijgsman et al. 2019). The following sea level rise is estimated for about 7,500 years BP (Cordova 2007; Fouache et al. 2012), which is consistent with the beginning of the Atlantic period. It is supposed that in this time, J. excelsa was abundant on the Crimean coastal zone (Yena et al. 2005), which could explain the ancestral character of its southern coastal populations in Crimea, proposed in our study based on the morphological similarities of these populations and the presence of the cpDNA H2 haplotype. A significant part of the Crimean coast with Mediterranean flora was eventually flooded, and the possibilities for plant migration around the Black Sea decreased (Yena et al. 2005). This scenario can be confirmed by the results of our phylogenetic analyses, which implied the common origination and recent fragmentation of populations of the species. These warm conditions could also foster the expansion of J. excelsa in Crimea northward and to higher locations in the mountains, which is consistent with the recent expansion model revealed by the minimum spanning network.