Skip to main content

Divergence of alpine plant populations of three Gentianaceae species in the Qinling sky Island

Abstract

Background

Known for their unique biodiversity, the Qinling Mountains are considered the only area in which alpine biomes occur in central China. Given that the alpine biomes are particularly sensitive to global warming, understanding how alpine plants respond to climatic fluctuations is essential for the evolution and conservation of biodiversity. To address this issue, three alpine species of the Gentianaceae (Gentiana crassuloides, G. hexaphylla and Swertia bifolia) that represent different life types and diverse genera were selected.

Results

Genetic clustering analysis according to around 33,317 to 185,133 SNPs showed that the Qinling population was a separate lineage within each species. A high level of genetic differentiation was observed among the Qinling populations and the other populations of each species. Divergence time estimation based on plastomes and approximate Bayesian computation based on genomic SNPs showed that Qinling populations of the three Gentianaceae species originated at different periods under various patterns including primary source and hybridization. Significant signals of isolation by distance and isolation by environment were found in all three species. The redundancy and gradient forest analyses revealed that several temperature- and precipitation-related variables mainly contributed to shaping the genetic differentiation among the Qinling populations and others, indicating that the three species exhibited a similar pattern of adaptations to local environments.

Conclusions

This study unveiled the unique genetic and evolutionary features of the Qinling populations of these three species and elucidated the contributing role of both the environmental gradient and geographical isolation in genetic differentiation, which scientifically supports future conservation efforts.

Peer Review reports

Introduction

Global large mountain systems, such as the region of the Qinghai-Tibet Plateau (QTP) and the Andes, are home to about one-third of the world’s terrestrial biodiversity and the key centers of endemism and species richness [1,2,, 3]. Many of these areas are qualified as hotspots of biodiversity that are under threat [45]. Expanding our knowledge of the responses of the remarkable biodiversity of assembled plant communities to climate change is a matter of utmost importance to evolutionary and conservation biology, in particular in the region known for its alpine flora, which is quite susceptible to global warming that drives habitat contraction and genetic erosion of alpine biodiversity [6,7,, 8].

The geological processes, e.g., uplift and erosion, and climatic variations, e.g., climate oscillations, play roles in the accumulation of mountain biodiversity over evolutionary time [1, 9, 10, 11]. For example, both the uplift and climatic changes caused the assembly of biodiverse plant communities in the QTP, which occurred via the local species recruitment and colonization events from adjacent regions [12]. Furthermore, diversification may take place in adaptation to climate-driven shifts in the distribution ranges of species. The elevational shift of the mountainous plant species and the range fragmentation and isolation in sky islands resulting from global warming consequently contribute to allopatric diversification [1314], and the secondary contact after allopatric divergence caused by the downslope movement of plants during climate cooling could enable the promotion of hybridization [15]. Meanwhile, montane species experience upslope shifts in the face of global warming, upon which habitat reduction occurs in most mountain systems [1617]. Therefore, a pressing concern is the exploration of the origin of alpine plants with a very limited capacity for upslope migration under continued global warming, in particular in the sky island, in order to contribute to their conservation [14]. In addition, as the ideal natural laboratory for exploring the evolutionary consequences of climate change, the comparison of sky islands with the adjacent mountain systems could provide significant insights into the driving mechanisms for ecological and evolutionary responses to a changing climate [14].

The Qinling Mountains (QM) is a typical sky island for alpine taxa, and its main peak, Mount Taibai, is the highest mountain in Central and Eastern China with an altitude of 3500 m a.s.l. It stretches more than 400 km on the northeast margin of the QTP, and is considered a widely recognized geological and geographical boundary between northern and southern China. Additionally, the QM is a major biodiversity hotspot in China [18] and is renowned for high biodiversity and endemism richness, e.g., animal species such as giant pandas, snub-nosed monkeys, and crested ibises. In fact, 5.6% (192) and 45.7% (1428) of plant species in the QM are endemic to the QM and China, respectively (Ying, 1994). The Qinling flora is a key part of the East Asian flora [19], but also contains rich alpine taxa. The alpine flora in the QM comprises not only the typical alpine taxa composed of widely distributed taxa such as fir, Gentiana, and Saxifraga but also the endemic species such as Oxytropis chinglingensis C.W. Chang, Petrocosmea qinlingensis W. T. Wang, and Larix potaninii var. chinensis L.K. Fu & Nan Li, pointing out local diversification and the colonization of the QM and its adjacent regions by alpine flora. The QTP, a potential target for colonization, is characterized by rich alpine biodiversity and the oldest alpine lineages [12]. Although the colonization of the QTP and the QM by alpine lineages is seemingly frequent, its associated factors, e.g., the direction of colonization and similar potential patterns of diverse taxa, have not yet been estimated well.

The QM region accommodates a typical alpine flora like Gentianaceae, a family composed of several species-rich alpine lineages such as Gentiana and Swertia, both of which were found to have the highest species diversity in the QTP region [20,21,, 22]. A total of 15 Gentianaceae species exist in the QM, among which 10 occur above 3000 m [23] (Table S1). Apart from one species endemic to the QM region (Gentiana apiata N.E.Br.), eight species distributed in both the QM and the QTP regions, such as Gentiana hexaphylla Maxim and Swertia bifolia Batalin [20, 23]. Therefore, Gentianaceae represents an ideal model taxon for exploring the evolutionary relationship between alpine plant species in both the QM and the QTP regions.

To better understand the underlying reasons and evolutionary processes of alpine species in the QM, three alpine species of Gentianaceae, including one annual and one perennial gentian species, as well as one perennial Swertia species, that are representatives of different life types and diverse genera, were selected. Although the populations of the three species had relative large population sizes in the QM, they undoubtedly would face the risk of habitat reduction during upslope shifts under global warming. The study aimed to address the following questions: (1) Do the QM populations of the Gentianaceae species possess unique genetic and evolutionary features? (2) Do Gentianaceae populations in the QM contemporarily diverge and act as the primary source? and (3) What are the causative factors for the genetic divergence between plant populations in the QM and adjacent regions? By combining genomic and bioclimatic data, our analyses could provide evidences on the genetic traits and origin of alpine plants in the QM, and shed new light on the biodiversity conservation in the Qinling sky island.

Results

Data preprocessing and SNP calling

There were variations in the number of retained clean reads per sample for G. crassuloides, ranging from 0.66 × 106 to 5.40 × 106, with an average value of 2.11 × 106, while for G. hexaphylla, it varied from 4.44 × 106 to 80.68 × 106, averaging at 16.78 × 106, and, the variation within the range of 1.97 × 106–5.11 × 106 and an average value of 3.13 × 106 were obtained for S. bifolia. The pruning of SNPs for minor allele frequency, missing data, and linkage disequilibrium resulted in the total number of unlinked SNPs of 33,569, 185,133, and 33,317 for G. crassuloides, G. hexaphylla, and S. bifolia, respectively, for all download analyses except DIYABC-RF. During DIYABC-RF analyses, since we only kept SNPs with no missing among individuals, 5,961, 37,273, 7,483 SNPs were remained for G. crassuloides, G. hexaphylla, and S. bifolia, respectively.

Population genetic structure and divergence

The analysis of genome-wide SNPs revealed a clear genetic structure of the three Gentianaceae species and a deep genetic divergence within each species. For G. crassuloides, the gTB population of the Qinling had a slightly lower value of HO (0.103) but a significantly higher value of FIS (0.347) than the QTP populations (average values = 0.112 and 0.272; p = 0.152 and 0.002, respectively) (Table 1). The fastSTRUCTURE analysis that maximized the model complexity and marginal likelihood score was K = 6 in G. crassuloides (Fig. S4). The six clusters were corresponded to their geological distribution (Fig. 1a and d). The gTB population of the Qinling formed a separate cluster. Similar to Bayesian clustering results, PCA plots also showed the grouping of all individuals into six clusters. PCA results demonstrated that the first two principal components (PC1 and PC2) explained 38.66% and 13.24% of the total genetic variance, respectively, indicating the complete differentiation of these clusters (Fig. 1b). Genetic distances (pairwise FST values) among clusters were much higher than those within clusters (Fig. 1c). The FST values between the gTB population and others varied from 0.431 to 0.714 (Table S3).

Table 1 Sample information and genetic statistics for the three Gentianaceae species
Fig. 1
figure 1

Geographic distributions and genetic differentiation of three Gentianaceae species based on nuclear SNP data. Colors indicate genetic groups clustered based on the fastSTRUCTURE software tool and principal coordinate analysis. (1) a-d: results of the population genetic analysis using 33,569 SNPs for Gentiana crassuloides. (a) Geographic distributions of the six clusters; (b) results of the principal coordinate analysis; (c) the heatmap of weighted Weir and Cockerham’s FST values; (d) the bar plots showing the probabilities of ancestral clusters of each sample at K = 6. (2) e-h: results of the population genetic analysis based on 185,133 SNPs for Gentiana hexaphylla. (e) Geographic distributions of the four clusters; (f) results of the principal coordinate analysis; (g) the heatmap of weighted Weir and Cockerham’s FST; (h) the bar plots display the assignment probabilities of individuals to distinct ancestral clusters for each sample at K = 4. (3) i-m: results of the population genetic analysis using 33,317 SNPs for Swertia bifolia. (i) Geographic distributions of the two clusters; (j) results of the principal coordinate analysis; (k) the heatmap showing pairwise weighted Weir and Cockerham’s FST values; (m) bar plots represent the probability of the assignment of each individual to each ancestral cluster at K = 2

A significantly lower HO value (0.138) but a significantly higher FIS value (0.241) were found for the TB population of G. hexaphylla from the Qinling than those of the QTP populations (average = 0.184 and 0.061; p = 3.83E− 7 and 3.89E− 7, respectively) (Table 1). Four genetic clusters for G. hexaphylla were detected by the Bayesian clustering method since the model complexity and marginal likelihood score was maximized at K = 4 (Fig. S4). The individuals from the TB population from the Qinling formed an independent cluster, and all clustering patterns corresponded to the geological distribution of the population (Fig. 1e and h). The cluster analysis based on PCA showed four well-differentiated clusters, with PC1 and PC2 explaining 8.59% and 5.02% of the total genetic variance, respectively (Fig. 1f). There were high genetic distances (FST values) between populations from different clusters (Fig. 1g). The FST values between the TB population and the others were within the range of 0.216–0.305 (Table S4).

With regards to S. bifolia, the sTB population of the Qinling had slightly lower HO (0.169) and FIS (0.148) values than the QTP populations (average values = 0.176 and 0.157; p = 0.377 and 0.789, respectively) (Table 1). According to Bayesian clustering, individuals of S. bifolia were divided into two clusters (K = 2; Fig. S4). The sTB population sampled from the QM formed a separate cluster, while all other populations belonged to the other (Fig. 1i and m). PCA could fully differentiate the two clusters, and 4.55% and 3.13% of the total genetic variance were explained by PC1 and PC2, respectively (Fig. 1j). The sTB population exhibited much higher FST values than the other populations (Fig. 1k), and there was a variation of 0.210 to 0.272 in FST between the sTB population and the others (Table S5).

Phylogenetic relationships and divergence time

Our genomic SNP data generated well-supported ML trees of the three Gentianaceae species, with 100% bootstrap support for most of the deep nodes (Fig. 2). The samples of G. crassuloides were divided into three main sub-lineages as revealed by the analysis of genetic clustering using fastSTRUCTURE and PCA, and the 4th sub-lineage was formed with the gTB population from the QM and found to be at the base of the tree. For G. hexaphylla, the TB population from the QM was sister to a clade containing the North and partial Central populations, which were themselves sister to a clade composed of the South populations. One Central population of G. hexaphylla was observed to be located at the base of the tree. The ML tree of S. bifolia showed that the sTB population from the QM located at the base of the tree and was sister to a sub-lineage involving samples from the West cluster.

Fig. 2
figure 2

Maximum likelihood trees of Gentiana crassuloides (a), G. hexaphylla (b), and Swertia bifolia (c) constructed based on genome-wide SNPs. The support values for maximum likelihood on phylogenetic trees are indicated near the main nodes, with asterisks showing 100% support. The colors correspond to the genetic clusters based on genome-wide SNPs

Our divergence time analyses based on plastomes showed that Gentianinae diverged from Swertiinae 60.84 Ma (95% HPD: 59.91–61.75 Ma). Divergence among the three sublineages of G. crassuloides occurred 3.13 Ma (95% HPD: 2.93–3.34 Ma) and 2.65 Ma (95% HPD: 2.46–2.83 Ma), respectively. The Qinling population gTB diverged from its sister clade 1.18 Ma (95% HPD: 1.05–1.32 Ma). In G. hexaphylla, the two sublineages diverged 1.72 Ma (95% HPD: 1.53–1.90 Ma). The Qinling population TB diverged from its sister clade 0.11 Ma (95% HPD: 0.06–0.17 Ma). Divergence within S. bifolia began at 0.79 Ma (95% HPD: 0.66–0.93 Ma), and the Qinling population sTB subsequently diverged from its sister clade 0.36 Ma (95% HPD: 0.28–0.45 Ma) (Fig. 3).

Fig. 3
figure 3

Divergence time estimation based on plastomes in the three Gentianaceae species. The gray bars show the 95% highest posterior density on the age estimates. The red circles show the two calibration nodes. The posterior probabilities of nodes were presented above branches. Ma, million years ago; QU, Quaternary; PL, Pliocene. The populations from the Qinling Mountains are marked in red

Colonization scenarios of the three Gentianaceae species

For all DIYABC-RF analyses, the random forest LDA projections of simulated data and error metrics fit the observed data, indicating that the model choices were reliable (Fig. S5). Our DIYABC-RF analysis in G. crassuloides supported model 11 with the largest number of votes and the highest posterior probability (Fig. S1, Table S6). The most likely scenario suggested an ancestral origin of the gTB population of the Qinling, followed by an early colonization of the northern and southern QTP, and subsequent dispersal to the western QTP through hybridization (Fig. 4a). Our optimal model estimated the occurrence of the early dispersal to the Qinling Mts. 1.35E6 (95% HPD: 410,961–2.07E6) generations (years) ago, and to the western QTP areas 280,397 (95% HPD: 161,419–628,855) generations (years) ago. The effective population size of the Qinling population was estimated to be 581,972 (Table S7). These demographic and temporal estimates for the gTB population were in general agreement with phylogenetic reconstruction (Fig. 2a) and dating analysis (Fig. 3).

Fig. 4
figure 4

Graphical representation and time parameter estimation of the selected scenarios on the alternative dispersal models of the population groups of the three Gentianaceae species from DIYABC-RF analysis. The right scale shows the states for the time parameter that generation time multiplies the median estimation. Migration rates were presented when hybridization events occurred. Colors showed the different population groups as the follow. (a) Gentiana crassuloides: North (gLX, gGZ, gMEK, gHY, gJZ1), West (gREG, gLH, gYS, gLWQ, gJZ2, gDR), South (gYJ, gLD, gDQ, gCY), and gTB. (b) Gentiana hexaphylla: North (JZ, HY), Center (SE, LH, KD), South (XC, DQ, GS, CY), and TB. (c) Swertia bifolia: QTP (sSD, sXH, sDR, sGD, sBM, sCD1, sCD2, sJD, sDQ), and sTB. Post. prob.: posterior probability

Regarding to G. hexaphylla, the DIYABC-RF analysis supported model 8 with the largest number of votes and the highest posterior probability (Fig. S1, Table S6). The most likely scenario suggested an ancestral origin of the populations of the southern QTP, followed by a colonization of the Qinling Mts. and the northern QTP, and subsequent dispersal to the central QTP through hybridization between southern and northern populations of the QTP (Fig. 4b). Setting the generation time as two years for the perennial herb, our optimal model estimated the occurrence of the early dispersal to the southern QTP about 2.36 Ma (95% HPD: 1.51–5.76 Ma), and to the Qinling Mts. 1.25 Ma (95% HPD: 0.22–5.16 Ma). The effective population size of the Qinling population was predicted to be 534,348 (Table S7).

In turn to S. bifolia, model 5 was the most likely scenario (Fig. S1, Table S6) and suggested an ancestral origin of the populations of the QTP, the sTB population originated through hybridization between the QTP populations and a ghost population (Fig. 4c). Setting the generation time as two years, the optimal model estimated the occurrence of the early dispersal to the QTP about 5.24 Ma (95% HPD: 0.51–19.62 Ma), and to the Qinling Mts. 0.25 Ma (95% HPD: 0.05–0.52 Ma). The effective population size of the Qinling population was estimated to be 695,373 (Table S7).

Environmental and spatial associations with the genetic differentiation

The results of the redundancy analysis (RDA) showed that significant proportions of genetic variation among the individuals of the three Gentianaceae species were shaped by bioclimatic variables, and the important predictors of genetic differentiation were different in the three species. For G. crassuloides, eight unrelated bioclimatic variables contributed to the great proportion of genetic variation (33.43%, p = 0.001), with all of these variables found to have significant associations with genetic differentiation (p = 0.001; Table 2). Annual precipitation (bio12, 9.48%), temperature annual range (bio7, 9.44%), and precipitation seasonality (bio15, 8.44%) were the three most important predictors of genetic variation (Table 2). The main variable explaining the differentiation between the East cluster composed of the gTB population and the other clusters was the annual mean temperature (bio1) (Fig. 5a). RDA indicated that nine bioclimatic variables had a significant influence (p = 0.001) on 28.12% of the genetic variation in G. hexaphylla. All nine variables were found to be significantly related to genetic variation (p = 0.001; Table 2), among which the precipitation of the wettest month (bio13, 9.48%), temperature seasonality (bio4, 6.21%), and temperature annual range (bio7, 5.73%) were the top three predictor variables (Table 2). Annual mean temperature (bio1), annual precipitation (bio12), and precipitation of the driest month (bio14) were the main factors causing the genetic differentiation between the TB population from the Qinling and the others (Fig. 5b). Concerning S. bifolia, 29.1% of the genetic variation was shown to be significantly affected by seven bioclimatic variables (p = 0.001) according to RDA analysis, all of which were found to have significant relationships with genetic variation (P = 0.001 or 0.003; Table 2). Among them, the annual precipitation (bio12, 5.33%), mean diurnal range (bio2, 5.04%) and max temperature of warmest month (bio5, 4.86%) were the three most important factors (Table 2). The differentiation between the sTB population from the Qinling and the others was mainly explained by the minimum temperature of the coldest month (bio6) and annual precipitation (bio12) (Fig. 5c).

Table 2 Results of the redundancy analysis (RDA) in the three Gentianaceae species based on un-related environmental variables
Fig. 5
figure 5

The redundancy analysis showing the relationship between unrelated climatic variables and the genetic divergence among populations of (a) Gentiana crassuloides, (b) G. hexaphylla, and (c) Swertia bifolia. The colored points represent the individuals of the genetically divergent clusters inferred from a nuclear SNP dataset. Bio1, annual mean temperature; bio2, mean diurnal range; bio3, isothermality; bio4, temperature seasonality; bio5, maximum temperature of the warmest month; bio6, minimum temperature of the coldest month; bio7, temperature annual range; bio12, annual precipitation; bio14, precipitation of the driest month; and bio15, precipitation seasonality; ***P < 0.001, ** P < 0.005

Our gradient forest (GF) analysis revealed that the key variables predicting allele frequency differed among the three Gentianaceae species, indicating strong turnover among the Qinling populations and the others (Fig. 6c). Annual precipitation (bio12) and temperature seasonality (bio4) were the top two variables that could predict allele frequency for both G. crassuloides and S. bifolia, while the two most important variables predicting allele frequency in G. hexaphylla were the temperature seasonality (bio4) and temperature annual range (bio7) (Fig. 6a). The abrupt change in the cumulative allele frequencies in G. crassuloides was recorded with the annual precipitation of 600–620 mm, temperature seasonality of 630–680, precipitation seasonality of 87–90, and annual mean temperature of -1–0 ℃. The cumulative allele frequencies in S. bifolia evolved across the environmental gradient and they changed suddenly with the annual precipitation of 630–660 mm and temperature seasonality of 750–770. However, at a temperature seasonality of 610–650, a temperature annual range of 28–29 ℃, and a mean diurnal range of about 14 ℃, abrupt changes in the cumulative allele frequencies in G. hexaphylla occurred (Fig. 6b).

Fig. 6
figure 6

Results of the gradient forest (GF) analysis demonstrate the patterns of genetic variation in the three Gentianaceae species. (a) The R2-weighted importance of environmental variables explaining the genetic gradients based on the GF analysis. (b) Cumulative importance of the change in allele frequency across the top four environmental gradients; bio1, annual mean temperature (℃); bio2, mean diurnal range; bio3, isothermality (×100); bio4, temperature seasonality (standard deviation ×100); bio5, maximum temperature of the warmest month; bio7, temperature annual range; bio12, annual precipitation (mm); and bio15, precipitation seasonality (the coefficient of variation). (c) The spatial turnover of allele frequencies predicted with all SNPs set as the response variables and environmental variables set as explanatory variables. The genetic turnover of each predicted variable is reduced into three principal components that was each assigned an RGB color. The similarity between colors signifies the GF-predictable adaptive genetic frequencies and their differences between populations. The arrow indicates the populations sampled from the Qinling Mountains

There was a high correlation between pairwise genetic distance and geographical distance in the three Gentianaceae species, suggesting strong IBD (p = 0.001 or 0.003; Fig. 7a). A significant correlation was also identified between the genetic distance and environmental distance in all these three species (p = 0.001 or 0.055; Fig. 7b). In addition, a strong autocorrelation was also shown between the environmental and geographical distances (p = 0.001 or 0.003; Fig. 7c). Based on the partial Mantel test, the significant impact of geographical distance in G. crassuloides (p = 0.01), G. hexaphylla (p = 0.001) and S. bifolia (p = 0.046) was noted, while environmental distance exerted no significant effect in the three species (Table S8), suggesting the stronger role of geographical distance in structuring intraspecific diversity.

Fig. 7
figure 7

Mantel tests of the genetic distance [FST/ (1- FST)] against the geographic distance (a) and the environmental distance (b), the correlation between the geographic and environmental distances (c) in the three Gentianaceae species. The Mantel’s r values (r) and the statistical significance (p-value) represented at the right bottom of each panel

Discussion

The QM is a top biodiversity hotspot in China and a typical sky island in the world that inhabits alpine flora whose suitable habitat is declining under climate change scenarios. Therefore, of critical interest is to explore the genetic traits and origin of alpine plants in the Qinling Mountains that underlie their conservation. Taking the three alpine Gentianaceae species as the research object, in this study, separate clusters were found to be formed with the Qinling populations, between which and the QTP populations, a deep genetic differentiation occurred. Both geologic isolation and climatic adaptation contributed to differentiation, and the relatively low genetic diversity in the Qinling populations resulted from both the high levels of inbreeding and low gene flow. Here, our findings were addressed from the standpoint of evolutionary biology, and the challenges of the biodiversity conservation of the Qinling populations were contemplated.

Deep genetic divergence between the Qinling and the QTP populations

As the region that harbors the world’s richest temperate alpine flora, the alpine biotas in the QTP were assembled by both the local recruitment and the colonization of adjacent regions [12]. Despite its rich alpine flora, there is an inadequate investigation of the genetic relationships between the populations of this region and Qinling populations of alpine plants. To compensate for this lack of discussion, the three Gentianaceae species that occur in both the QTP and the QM regions were selected, and a large number of genome-wide SNPs were applied in this study. According to our results, the Qinling populations of all three species were grouped into a separate cluster, and there was a great genetic differentiation between these populations and the QTP populations, with FST values ranging from 0.236 to 0.539. Similarly, genetic divergence existed between the Qinling populations of some alpine plants such as Allium cyaneum [24], Notopterygium incisum [25], Pinus armandii [26], and Rheum palmatum [27] and the QTP populations. The isolation of sky islands also generated the lineages of floras in the Andes [28]. However, the QM seems to have a more complex situation, e.g., the Qinling populations of the two alpine herbs, Circaeaster agrestis [29] and Allium sikkimense [24], formed a genetic cluster with populations from the eastern QTP. The diverging lineages of an alpine herb, Notopterygium oviforme [30], and a woody plant in a montane forest Dipteronia sinensis [31] have even been present in the QM. Thus, the genetic relatedness in the Qinling and QTP populations of plants should supposedly be species-specific, but there was an inter-population genetic differentiation between alpine plants, especially the Gentianaceae family, from the QM and the QTP.

Besides an evident genetic divergence between the Qinling and the QTP populations, a distinct genetic differentiation was also detected among the populations of the three Gentianaceae species from the QTP. Great genetic distances (high FST values) observed among the clusters were identified by both the genetic clustering and PCA. Indeed, a distinct genetic differentiation within a species occurring in high mountain ranges is generally associated with geological isolation and high environmental heterogeneity [29, 32]. In addition, a stronger intraspecific differentiation of the annual gentian species G. crassuloides (FST = 0.475) than that of the perennial species G. hexaphylla (FST = 0.211) and S. bifolia (FST = 0.158) was found. This is in line with the hypothesis that annual herbs tend to have a higher level of genetic divergence as well as higher rates of molecular evolution than perennials [33] and that the role of life history differences in predicting genomic divergence. Moreover, in his study, the Qinling populations of the two Gentiana species, in particular, G. hexaphylla, an outcrossing dominant species [34] were found to have a significantly lower genetic diversity or/and a higher inbreeding coefficient than the QTP populations. The geological barrier in a sky island resulted in both a high inbreeding coefficient and low gene flow, which together were the contributory factors of the relatively low genetic diversity in Qinling populations of the two gentians. The reduced genetic variation among the populations of the sky island and the populations of the adjacent regions was noted [14].

The contribution of climatic and geographic factors to genetic differentiation

The geographical isolation, geological factors, climatic fluctuations, and local adaptation were the factors determining the intraspecific genetic differentiation [35,36,37,, 38] and varied among taxa [3940]. Our correlation analysis showed a significant correlation between the genetic and geographic distances in all three Gentianaceae species. Although the genetic and environmental distance were also significantly correlated, the partial Mantel test revealed the key role of geographic distance. These results indicated that IBD as well as IBE were two factors associated with the spatial patterns of intraspecific divergence in the three species in the QM and QTP, consistent with the findings obtained in mountain ranges including the Alps [41] and the Andes [4243]. Furthermore, the redundancy analysis confirmed the pivotal roles of climatic factors, including temperature-related variables, e.g., annual mean temperature (bio1), and precipitation-related variables, e.g., annual precipitation (bio12), in shaping genetic differentiation between the Qinling and the QTP populations. There was also a strong genetic turnover in the Qinling and QTP populations, and the major variables predicting allele frequency were the precipitation-related factor, e.g., bio12, and the temperature-related factor, e.g., temperature seasonality (bio4) as revealed by the GF analysis. These results revealed the climatic adaptation of the three Gentianaceae species, a case similar to that of montane biota in the Alps [4445] and the Andes [46]. The reason for the differences in the major climatic factors that caused genetic differentiation between the Qinling and QTP populations in the three species may be the difference in the preferred habitat, for example scrubs for S. bifolia and wet habitat including stream banks and bogs for G. crassuloides [23]. Altogether, the geographical and environmental factors were seen to be of importance in shaping the intraspecific diversification in the three Gentianaceae species, with a comparable finding reported for the Qinling species, including Circaeaster agrestis [29], Notopterygium oviforme [30], and Pterocarya macroptera [47]. Taking this into account together with the alpine endemic species such as Gentiana apiata and Oxytropis chinglingensis, these results indicated that the diversification caused by geographical isolation and local adaptation in the QM, e.g., the Taibai Mountain, as a sky island famous for its alpine flora, should be accelerated.

The diverse origin of the populations of gentianaceae species in the Qinling region

Although the alpine flora of the Qinling and the QTP regions have supposedly close relationships, their origin in the former remains to be clarified. Results from previous studies reported that the origin of the Qinling populations varied among different plants. For example, the Qinling population of species such as Paeonia qiui [48], Pinus armandii [26], and Rheum palmatum [27] and populations from the Daba Shan area shared a similar genetic composition, and likewise, there was a genetic link between the population of Circaeaster agrestis [29] in the Qinling region and that from the eastern QTP. A study conducted by Wang et al. [47] suggests a potential hybrid origin of Qinling populations of the woody plant in a montane forest ecosystem, Pterocarya macroptera.

Taking the three alpine herbs as the focus of this study, our new findings indicate that the QM could be a primary origin of alpine plants, in particular G. crassuloides, supported by both phylogenetic relationship (Fig. 2) and evolutionary history inferred by ABC analysis (Fig. 4a). Divergence time estimation based on plastome and ABC analysis based on nuclear SNPs consistently indicated that the gTB Qinling population diverged at the early Pleistocene when significant climate fluctuation including warm/cold and humid/arid alternating occurred in the Qinling Mts [49]. Interestingly, the TB Qinling population of G. hexaphylla also diverged from its sister group at the early Pleistocene (Fig. 4b). These evidences from the two species consistently showed that the climate fluctuation shall be one key driver for the genetic differentiation for the two gentians. Additionally, the gentian species G. apiata endemic to the QM should have arisen by geological isolation, highlighting the in situ differentiation of Gentiana species in this region. Results of S. bifolia showed that the sTB Qinling population originated from hybridization, indicates the diverse origin of Gentianaceae species in the QM. Hybridization events revealed by ABC analysis (Fig. 4) and phylogenetic discordance between nuclear and chloroplast genomes (Figs. 2 and 3) clearly showed that hybridization shall be common in the target species. Undoubtly, growing evidences showed that hybridization may be rampant in Gentiana [32, 50, 51, 52, 53] and alpine plants in the QTP [54]. In summary, the current limited evidence suggests that the origin of the alpine populations of the Qinling Mountains varied based on taxa and that the QM is an important region where partial alpine plants occurring in both the QM and QTP come from.

Implications for biodiversity conservation

Identifying different genetic units by the determination of the boundary of the unit for policymakers is a crucial step in the conservation of biodiversity. The unique genetic composition and evolutionary features of the Qinling populations of the three Gentianaceae species were shown in the present study. The QM particularly acted as the original hotspot of alpine plants such as G. crassuloides and S. bifolia. We observed that the Qinling populations shall not be at risk in the wild due to current large population sizes, and our estimations indeed showed large effective population sizes. However, the survival of the three alpine Gentianaceae species in the QM is under threat of habitat contraction induced by climate warming [16] and human activities such as continuous harvesting, the rise of hiking, and tourism activities. Although natural reserves have been created in the QM, the biodiversity of montane ecosystems comprised of animals such as giant pandas and snub-nosed monkeys seemingly receives more attention than that of the alpine flora. Therefore, a greater emphasis should be placed on designing scientific management practices and conservation measures for the protection of natural resources of the three Gentianaceae species. The protection of alpine habitats in the Qinling Mountains occupied by these species is a prerequisite to the restoration of wildlife resources. Indeed, protecting such habitats is beneficial to not only the three target species but also an endemic species, Gentiana apiata, and an endangered species Notopterygium oviforme [23] in alpine habitats of the QM. The implementation of ex situ conservation methods, such as transplanting into alpine botanical gardens, in addition to in situ conservation, should be taken into consideration in the near future. For G. hexaphylla and G. crassuloides, which have significantly high levels of inbreeding and low genetic diversity in the Qinling populations, a particularly essential process is transplanting. During ex situ conservation, besides the key climatic factors revealed in this study, soil characters such as hydrology in the primary regions [5556] shall also be considered.

Materials and methods

Study species and sampling

The present study focused on three alpine Gentianaceae species, including Gentiana crassuloides Franch., an annual herb, and Gentiana hexaphylla and Swertia bifolia, two perennial herbs, distributing in both the QM and the QTP regions. In details, G. hexaphylla and S. bifolia were endemic to China, distributing in Shaanxi, Gansu, Qinghai, Sichuan, Yunnan and Xizang; G. crassuloides overlapped its range in China with G. hexaphylla and S. bifolia, and also had range in Bhutan, NW India and Nepal [2021]. The presence of components such as gentiopicroside and swertiamarin contributes to the significant medicinal value of these three species [2021, 57]. A total of 89 and 44 individuals of G. crassuloides and S. bifolia, respectively, were freshly collected to show their ranges. The subsamples of 59 individuals of G. hexaphylla were retrieved from the study conducted by Fu et al. [58]. Altogether, 192 individuals in 10–16 different populations of the three species were sampled (Table 1). The voucher specimens of plant materials identified by Dr. Peng-Cheng Fu were deposited in the herbarium of Luoyang Normal University. The detailed information on specimens of the three species was presented in Table S2.

Illumina sequencing and data processing

The genotype-by-sequencing (GBS) method was employed to obtain single nucleotide polymorphism (SNP) genotypes for collected samples of G. crassuloides and S. bifolia, following the method described by Elshire et al. [59]. Total genomic DNA was extracted from dry leaves, and thereafter, the construction of libraries according to the procedure described in Fu et al. [32] were performed, followed by their sequencing using a 2 × 150 bp read length on two lanes of Illumina NovaSeq 6000 (Tianjin, China). To study population genetics in G. hexaphylla, data obtained by the restriction site-associated DNA sequencing (RAD-seq) method were retrieved from the research conducted by Fu et al. [58]. To remove adaptor sequences and low sequencing quality reads and bases, quality trimming of the raw reads was carried out using Trimmomatic v0.32 [60] with default parameters. The genome-alignment software bwa-men v2.2.1 [61] was used to map the clean reads against the available genome of the most closely related species, and the sequence alignment/map (SAM) format files were generated with SAMtools v1.9 [62]. We used the chromosome-level genome assemblies of Gentiana dahurica (GenBank no. PRJNA799480) [63] and Sinoswertia tetraptera (CNSA no. CNP0002339, https://db.cngb.org/cnsa) [64] as the reference genome for Gentiana and Swertia species samples, respectively. The PCR duplicate reads were marked using sambamba v0.8.1 [65], and variant calling was performed using freebayes v0.9.21 [66]. Raw variant calls were filtered using vcftools v0.1.13 [67], retaining only SNPs for the subsequent downstream analysis. This program was further used to remove the SNPs with minor allele frequency (MAF) values lower than 5%, missing frequency of higher than 0.8 among individuals, and size of less than 100 bp.

Genetic divergence and population genetic structure

The parameters such as allelic richness (Ar), observed heterozygosity (HO), gene diversity (HS) within populations, and Wright’s inbreeding coefficient (FIS) were computed using the R package Hierfstat [68]. The same package was also used to estimate the genetic differentiation among populations by the fixation index (FST) [69]. To visualize the results, the pheatmap package in R v.4.0.1 was used. We performed t-test to check if the differences of FIS and HO between the Qinling and QTP populations were significant in Excel.

Furthermore, we used the Bayesian clustering method implemented in fastSTRUCTURE [70] to investigate the genetic clusters in populations with the assumed number of clusters (K) ranging from 1 to 10. To assess model complexity for the data, chooseK.py script of fastSTRUCTURE software was employed, and DISTRUCT v1.1 software [71] was further used to visualize the resultant cluster assignments. The principal component analysis (PCA) was performed on the same data to investigate genetic structure among populations, and 10 main principal components (PCs) were obtained in PLINK v1.90 [72].

Phylogenetic reconstruction and divergence time estimation

The phylogenetic trees were built with the genome-wide SNPs of species based on the maximum likelihood (ML) method using the GRT nucleotide substitution model and 1000 ultrafast bootstrap replicates [73] implemented in IQ-TREE v.1.6.8 [74]. To convert the SNP data for use as the input for phylogenetic tree construction, we used a customized Python script ‘vcf2phylip’ (https://github.com/edgardomortiz/vcf2phylip). Gentiana aristata (Gao2019688), Gentiana hoae (Fu2017046), and Swertia wardii (Fu2016194) were the outgroup species selected according to previous phylogenetic studies [75,76,, 77] for G. crassuloides, G. hexaphylla, and S. bifolia, respectively.

The divergence times of the main lineages of target species were estimated using the Bayesian method implemented in BEAST 2.4 [7879]. The analyses were based on the plastome (removed one IR) matrix, and were performed using the GTR substitution, the calibrated Yule, and the strict clock model. As a result of poor fossil in Gentianaceae, the dating was calibrated in two ways to increase its accuracy following Fu et al. [77]. First, we constrained the stem node of G. section Cruciata with a robust Cruciata fossil from the early Miocene [80]. Lognormal priors with an offset at 16.0 Ma, a mean of 1, and a standard deviation of 1.0 were used. Second, we also constrained the crown age of Gentiana after including primary lineages of this genus. We used uniform priors [81] with a lower age of 21 Ma and an upper age of 38 Ma to integrate the entire 95% Highest Posterior Density (HPD) from Janssens et al. [82]. Three independent Monte Carlo Markov chains (MCMC) with 10 million generations, sampling every 1000th generation, were run with the initial 20% discarded as burn-in. Effective Sample Size (ESS) values (> 200) were utilized to determine whether the convergence was appropriate, and TreeAnnotator 1.7.5 [78] was used to summarize the tree.

Population evolutionary relationships and demographic history

We evaluated the support for different possible scenarios of population divergence and admixture using t coalescence-based Approximate Bayesian Computation (ABC) and supervised machine learning methods implemented in DIYABC-RF (Random Forest) v.1.1.1 [83]. To reduce the computing load, we kept SNPs with no missing among individuals (--max-missing = 1) using vcftools v0.1.13 [67]. To trace the origin of the Qinling populations, we statistically tested alternative scenarios for the dispersal and admixture of populations of the three Gentianaceae species. We defined various main genomic groups in the three species identified by the fastSTRUCTURE analysis and according to geographical range. In G. crassuloides, we defined four main genomic groups [(i) West: gREG, gLH, gYS, gLWQ, gJZ2, gDR; (ii) North: gLX, gGZ, gMEK, gHY, gJZ1; (iii) South: gYJ, gLD, gDQ, gCY; (iv) TB: gTB] and assessed 12 competing scenarios for the dispersal history of the G. crassuloides groups based on our previous population genomic results, including the four groups originated from one ancient population one after the another (scenario 3–9), from two independent ancient populations (scenario 1–2), and one group originated from hybridization (scenario 10–12) (Fig. S1). In G. hexaphylla, we defined four main genomic groups [(i) North: JZ, HY; (ii) Center: SD, LH, KD; (iii) South: XC, CY, DQ, GS; (iv) TB: TB] and assessed eight competing scenarios based on our previous population genomic results, including the four groups originated from one ancient population at the same time (scenario 1) or one after another (scenario 5–7), from two independent ancient populations (scenario 2–4), and one group originated from hybridization (scenario 8) (Fig. S2). Regarding to S. bifolia, we defined two main genomic groups [(i) QTP: sSD, sXH, sDR, sGD, sBM, sCD1, sCD2, sJD, sDQ; (ii) TB: sTB] and assessed five scenarios for the S. bifolia groups, including the two groups originated from one ancient population (scenario 1), TB originated from QTP (scenario 2), QTP originated from TB (scenario 3), TB originated from hybridization between QTP and a ghost group (scenario 4), QTP originated from hybridization between TB and a ghost group (scenario 5) (Fig. S3). For all scenarios for each species, training sets were generated using 10,000 simulations per model. Due to the lack of ancestral divergence time information, we set the prior distribution of tree divergence times between the oldest and most recent splits of the coalescent-dated tree and the effective population sizes equal for each scenario of the same species. The most likely scenario was determined using the RF module of DIYABC-RF, computing a forest of 5,000 random trees and using the median with a 90% credibility interval as a point estimate. The model with the highest classification vote and the highest posterior probability was considered the most suitable. We estimated parameter values by running 10,000 simulations for the final selected model to infer time parameters and effective population size.

Environmental variables driving genetic divergence

To determine the relative contributions of environmental variables to shaping genetic differentiation in the three Gentianaceae species, we performed the gradient forest (GF) and redundancy analyses (RDA). GF is a machine learning algorithm that allows the prediction of non-linear relationships among the input spatial, environmental, and allelic variables that play a role in structuring the overall genetic variation [84]. A total of 19 bioclimatic variables were derived from the WorldClim database (https://www.worldclim.org/). One set of the variables that were found to be highly correlated (Pearson correlation coefficient > 0.8 and p < 0.01) using SPSS 20.0, were removed to avoid the occurrence of multicollinearity. The correlated climatic variables including 8 (bio1–bio4, bio7, bio12, bio14, and bio15), 9 (bio1–bio4, bio7, and bio12–bio15), and 7 (bio2–bio6, bio12, and bio14) were retained in the modeling of G. crassuloides, G. hexaphylla, and S. bifolia, respectively. We used the R package gradientForest v.0.1–34 to perform the GF analysis [85]. SNP turnover in allele frequencies was predicted by the GF model with all SNPs set as the response variables and environmental factors set as explanatory variables. The proportion of total genetic variance explained by all environmental variables was estimated by performing the RDA analysis with the function rda implemented in VEGAN v.2.5 [86] with the variables adopted by GF. To obtain the significance of the tested bioclimatic variables, we used the anova.cca function with 9999 permutations.

Inferring genetic differentiation patterns by isolation-by-distance and isolation-by-environment

We inferred the measures of genetic differentiation by isolation by spatial distance (IBD) and isolation by environment (IBE) to investigate how geographical and environmental factors influence spatial genetic differentiation. The geographical distance between pairs of populations was computed using the pointDistance function from the raster package [87] in R. The uncorrelated environmental variables were standardized to compare the relative magnitude of the effects of different explanatory variables and then used to calculate the environmental distance using the bcdist function from the ECODIST v.2.0.9 package in R [88]. The Mantel and partial Mantel tests [89] were performed to evaluate whether significant correlations existed among genetic distance (FST/(1-FST)), geographical distance (ln(km)), and environmental distance based on 999 permutations using R package VEGAN v.2.5 [86], whose plots were then created in R.

Data availability

The dataset supporting the conclusions of this article is available in the Figshare (https://doi.org/10.6084/m9.figshare.25671066). The raw sequence data reported in this paper have been deposited in the Genome Sequence Archive in National Genomics Data Center, China National Center for Bioinformation/Beijing Institute of Genomics, Chinese Academy of Sciences (GSA: CRA020255) that are publicly accessible at https://ngdc.cncb.ac.cn/gsa. The deposition number of plant samples was presented in Table S9.

References

  1. Antonelli A, Kissling WD, Flantua SGA, Bermúdez MA, Mulch A, Muellner-Riehl AN, et al. Geological and climatic influences on mountain biodiversity. Nat Geosci. 2018;11:718–25. https://doi.org/10.1038/s41561-018-0236-z.

    Article  CAS  Google Scholar 

  2. Noroozi J, Talebi A, Doostmohammadi M, Rumpf SB, Linder HP, Schneeweiss GM. Hotspots within a global biodiversity hotspot–areas of endemism are associated with high mountain ranges.Sci. Rep. 2018;8:10345. https://doi.org/10.1038/s41598-018-28504-9.

    Article  CAS  Google Scholar 

  3. Perrigo A, Hoorn C, Antonelli A. Why mountains matter for biodiversity. J Biogeogr. 2020;47:315–25. https://doi.org/10.1111/jbi.13731.

    Article  Google Scholar 

  4. Myers N, Mittermeier RA, Mittermeier CG, da Fonseca GAB, Kent J. Biodiversity hotspots for conservation priorities. Nature. 2000;403:853–8. https://doi.org/10.1038/35002501.

    Article  CAS  PubMed  Google Scholar 

  5. Bellard C, Leclerc C, Leroy B, Bakkenes M, Veloz S, Thuiller W, et al. Vulnerability of biodiversity hotspots to global change. Glob Ecol Biogeogr. 2014;23:1376–86. https://doi.org/10.1111/geb.12228.

    Article  Google Scholar 

  6. Pepin N, Bradley RS, Diaz HF, Baraer M. Elevation-dependent warming in mountain regions of the world. Nat Clim Change. 2015;5:424–30. https://doi.org/10.1038/nclimate2563.

    Article  Google Scholar 

  7. Meza-Joya FL, Morgan-Richards M, Koot EM, Trewick SA. Global warming leads to habitat loss and genetic erosion of alpine biodiversity. J Biogeogr. 2023;50:961–75. https://doi.org/10.1111/jbi.14590.

    Article  Google Scholar 

  8. Carruthers T, Moerland MS, Ebersbach J, Vaughan T, Wu CH, Xie D, et al. Repeated upslope biome shifts in Saxifraga during late-cenozoic climate cooling. Nat Commun. 2024;15:1100. https://doi.org/10.1038/s41467-024-45289-w.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Badgley C, Smiley TM, Terry R, Bermúdez MA, Mulch A, Muellner-Riehl AN, et al. Biodiversity and topographic complexity: modern and geohistorical perspectives. Trends Ecol Evol. 2017;32:211–26. https://doi.org/10.1016/j.tree.2016.12.010.

    Article  PubMed  PubMed Central  Google Scholar 

  10. Muellner-Riehl AN, Schnitzler J, Kissling WD, Mosbrugger V, Rijsdijk KF, Seijmonsbergen AC, et al. Origins of global mountain plant biodiversity: testing the ‘mountain-geobiodiversity hypothesis’. J Biogeogr. 2019;46:2826–38. https://doi.org/10.1111/jbi.13715.

    Article  Google Scholar 

  11. Li Q, Sun H, Boufford DE, Fennell T, Ruan J, Homer N, et al. Grade of membership models reveal geographical and environmental correlates of floristic structure in a temperate biodiversity hotspot. New Phytol. 2021;232:1424–35. https://doi.org/10.1111/nph.17443.

    Article  PubMed  Google Scholar 

  12. Ding WN, Ree RH, Spicer RA, Xing YW. Ancient orogenic and monsoon-driven assembly of the world’s richest temperate alpine flora. Science. 2020;369:578–81. https://doi.org/10.1126/science.abb4484.

    Article  CAS  PubMed  Google Scholar 

  13. He K, Jiang X. Sky islands of southwest China. I: an overview of phylogeographic patterns. Chin Sci Bull. 2014;59:585–97. https://doi.org/10.1007/s11434-013-0089-1.

    Article  Google Scholar 

  14. Love SJ, Schweitzer JA, Woolbright SA, Wang XQ. Sky Islands are a global tool for predicting the ecological and evolutionary consequences of climate change. Annu Rev Ecol Evol Syst. 2023;54:219–36. https://doi.org/10.1146/annurev-ecolsys-102221-050029.

    Article  Google Scholar 

  15. Ma YZ, Mao XX, Wang J, Zhang L, Jiang YZ, Geng YY, et al. Pervasive hybridization during evolutionary radiation of Rhododendron Subgenus hymenanthes in mountains of southwest China. Natl Sci Rev. 2022;9:nwac276. https://doi.org/10.1093/nsr/nwac276.

    Article  PubMed  PubMed Central  Google Scholar 

  16. Elsen PR, Tingley MW. Global mountain topography and the fate of montane species under climate change. Nat Clim Change. 2015;5:772–6. https://doi.org/10.1038/nclimate2656.

    Article  Google Scholar 

  17. Liang Q, Xu X, Mao K, Kang MH, Yang WJ, Feng LD, et al. Shifts in plant distributions in response to climate warming in a biodiversity hotspot, the Hengduan Mountains. J Biogeogr. 2018;45:1334–44. https://doi.org/10.1111/jbi.13229.

    Article  Google Scholar 

  18. Liu Q, Xue TT, Zhang XX, Yang XD, Qin F, Zhang WD, et al. Distribution and conservation of near threatened plants in China. Plant Divers. 2023b;45:272–83. https://doi.org/10.1016/j.pld.2023.02.005.

    Article  PubMed  PubMed Central  Google Scholar 

  19. Ying TS. An analysis of the flora of Qinling Mountain range: its nature, characteristics, and origins. J Syst Evol. 1994;32:389–410.

    Google Scholar 

  20. Ho TN, Liu SW. A Worldwide Monograph of Gentiana. Beijing: Science; 2001.

    Google Scholar 

  21. Ho TN, Liu SW. A worldwide monograph of Swertia and its allies. Beijing: Science; 2015.

    Google Scholar 

  22. Favre A, Michalak I, Chen CH, Albers CA, Banks E, DePristo MA, et al. Out-of-Tibet: the spatio-temporal evolution of Gentiana. (Gentianaceae) J Biogeogr. 2016;43:1967–78. https://doi.org/10.1111/jbi.12840.

    Article  Google Scholar 

  23. Ho TN, Pringle JS. Gentianaceae. In: Wu, Z.Y., Raven, P.H, editors. Flora of China, vol. 16. Beijing: Science Press; St. Louis: Missouri Botanical Garden. 1995; 1–140.

  24. Xie C, Xie DF, Zhong Y, Guo XL, Liu Q, Zhou SD, et al. The effect of Hengduan Mountains Region (HMR) uplift to environmental changes in the HMR and its eastern adjacent area: tracing the evolutionary history of Allium section Sikkimensia (Amaryllidaceae). Mol. Phylogenet. 2019;130:380–96. https://doi.org/10.1016/j.ympev.2018.09.011.

    Article  Google Scholar 

  25. Liu ML, He YL, López-Pujol J, Jia Y, Li ZH. Complex population evolutionary history of four cold-tolerant Notopterygium herb species in the Qinghai-Tibetan Plateau and adjacent areas. Heredity. 2019a;123:242–63. https://doi.org/10.1038/s41437-019-0186-2.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Liu YY, Jin WT, Wei XX, Wang XQ. Cryptic speciation in the Chinese white pine (Pinus armandii): implications for the high species diversity of conifers in the Hengduan Mountains, a global biodiversity hotspot. Mol Phylogenet. 2019;138:114–25. https://doi.org/10.1016/j.ympev.2019.05.

    Article  Google Scholar 

  27. Feng L, Ruhsam M, Wang YH, Li ZH, Wang XM. Using demographic model selection to untangle allopatric divergence and diversification mechanisms in the Rheum palmatum complex in the Eastern Asiatic Region. Mol Ecol. 2020;29:1791–805. https://doi.org/10.1111/mec.15448.

    Article  PubMed  Google Scholar 

  28. Vásquez DLA, Balslev H, Hansen MM, Sklenář P, Romoleroux K. Low genetic variation and high differentiation across sky island populations of Lupinus alopecuroides (Fabaceae) in the Northern Andes. Alp Bot. 2016;126:135–42. https://doi.org/10.1007/s00035-016-0165-7.

    Article  Google Scholar 

  29. Zhang X, Sun YX, Landis JB, Zhang JW, Yang LS, Lin N, et al. Genomic insights into adaptation to heterogeneous environments for the ancient relictual circaeaster agrestis (Circaeasteraceae, Ranunculales). New Phytol. 2020;228:285–301. https://doi.org/10.1111/nph.16669.

    Article  CAS  PubMed  Google Scholar 

  30. Liu ML, Shang QH, Cheng YJ, Jia Y, Li ZH, et al. Drivers of intraspecific differentiation of an alpine cold-tolerant herb, Notopterygium oviforme: roles of isolation by distance and ecological factors. J Syst Evol. 2023a;61:383–98. https://doi.org/10.1016/j.ympev.2019.05.015.

    Article  Google Scholar 

  31. Shahzad K, Liu ML, Zhao YH, Zhang TT, Liu JN, Li ZH. Evolutionary history of endangered and relict tree species Dipteronia sinensis in response to geological and climatic events in the Qinling Mountains and adjacent areas. Ecol Evol. 2020;10(24):14052–66. https://doi.org/10.1002/ece3.6996.

    Article  PubMed  PubMed Central  Google Scholar 

  32. Fu PC, Guo QQ, Chang D, Gao QB, Sun SS. Cryptic diversity and rampant hybridization in annual gentians on the Qinghai-Tibet Plateau revealed by population genomic analysis. Plant Divers. 2024;46:194–205. https://doi.org/10.1016/j.pld.2023.10.004.

    Article  PubMed  Google Scholar 

  33. Friedman J. The evolution of annual and perennial plant life histories: ecological correlates and genetic mechanisms. Annu Rev Ecol Evol Syst. 2020;51:461–81. https://doi.org/10.1146/annurev-ecolsys-110218-024638.

    Article  Google Scholar 

  34. He JD, Xue JY, Gao J, Wang JN, Wu Y. Adaptations of the floral characteristics and biomass allocation patterns of Gentiana hexaphylla to the altitudinal gradient of the eastern Qinghai-Tibet Plateau. J Mt Sci. 2017;14:1563–76. https://doi.org/10.1007/s11629-017-4424-x.

    Article  Google Scholar 

  35. Hewitt G. The genetic legacy of the Quaternary ice ages. Nature. 2000;405:907–13. https://doi.org/10.1038/35016000.

    Article  CAS  PubMed  Google Scholar 

  36. Rundle HD, Nosil P. Ecological speciation. Ecol Lett. 2005;8:336–52. https://doi.org/10.1111/j.1461-0248.2004.00715.x.

    Article  Google Scholar 

  37. Faurby S, Jørgensen A, Kristensen RM, Funch P. Distribution and speciation in marine intertidal tardigrades: testing the roles of climatic and geographical isolation. J Biogeogr. 2012;39:1596–607. https://doi.org/10.1111/j.1365-2699.2012.02720.x.

    Article  Google Scholar 

  38. Igea J, Tanentzap AJ. Global topographic uplift has elevated speciation in mammals and birds over the last 3 million years. Nat Ecol Evol. 2021;5:1530–5. https://doi.org/10.1038/s41559-021-01545-6.

    Article  PubMed  PubMed Central  Google Scholar 

  39. White OW, Reyes-Betancort JA, Chapman MA, Carine MA. Geographical isolation, habitat shifts and hybridisation in the diversification of the Macaronesian endemic genus Argyranthemum (Asteraceae). New Phytol. 2020;228:1953–71. https://doi.org/10.1111/nph.16980.

    Article  PubMed  Google Scholar 

  40. Hu HY, Yang YZ, Li A, Minh BQ, Vinh LS. Genomic divergence of stellera chamaejasme through local selection across the Qinghai–Tibet plateau and northern China. Mol Ecol. 2022;31:4782–96. https://doi.org/10.1111/mec.16622.

    Article  CAS  PubMed  Google Scholar 

  41. Wahlsteen E, Avramidou EV, Bozic G, Mediouni RM, Schuldt B, Sobolewska H. Continental-wide population genetics and post-pleistocene range expansion in field maple (Acer campestre L.), a subdominant temperate broadleaved tree species. Tree Genet Genomes. 2023;19(2):15. https://doi.org/10.1007/s11295-023-01590-1.

    Article  CAS  Google Scholar 

  42. Páez-Vacas MI, Trumbo DR, Funk WC. Contrasting environmental drivers of genetic and phenotypic divergence in an andean poison frog (Epipedobates anthonyi). Heredity. 2022;128:33–44. https://doi.org/10.1038/s41437-021-00481-2.

    Article  CAS  PubMed  Google Scholar 

  43. Muñoz-Valencia V, Montoya-Lerma J, Seppä P, Diaz F. Landscape genetics across the Andes mountains: environmental variation drives genetic divergence in the leaf-cutting ant Atta cephalotes. Mol Ecol. 2023;32:95–109. https://doi.org/10.1111/mec.16742.

    Article  PubMed  Google Scholar 

  44. Ceresa F, Brambilla M, Kvist L, Vitulano S, Pes M, Tomasi L, et al. Restricted dispersal and inbreeding in a high-elevation bird across the ‘sky islands’ of the European Alps. J Biogeogr. 2023;51:853–68. https://doi.org/10.1111/jbi.14787.

    Article  Google Scholar 

  45. Felkel S, Tremetsberger K, Moser D, Dohm JC, Himmelbauer H, Winkler M. Genome-environment associations along elevation gradients in two snowbed species of the North-Eastern Calcareous Alps. BMC Plant Biol. 2023;23:203. https://doi.org/10.1186/s12870-023-04187-x.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Sekely J, Marchelli P, Arana V, Rumpf SB, Linder HP, Schneeweiss GM, et al. Genomic responses to climate: understanding local adaptation in the Andean tree species Nothofagus Pumilio and implications for a changing world. Plants People Planet. 2024;6:902–20. https://doi.org/10.1002/ppp3.10504.

    Article  Google Scholar 

  47. Wang TR, Meng HH, Wang N, Zheng SS, Jiang Y, Lin DQ, et al. Adaptive divergence and genetic vulnerability of relict species under climate change: a case study of Pterocarya macroptera. Ann Bot. 2023;132:241–54. https://doi.org/10.1093/aob/mcad083.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Xu XX, Cheng FY, Peng LP, Sun YQ, Hu XG, Li SY, et al. Late pleistocene speciation of three closely related tree peonies endemic to the Qinling–Daba Mountains, a major glacial refugium in Central China. Ecol Evol. 2019;9(13):7528–48. https://doi.org/10.1002/ece3.5284.

    Article  PubMed  PubMed Central  Google Scholar 

  49. Lu HY, Zhang HY, Wang SJ, Zhao CF, Zhao J. A preliminary survey on loess deposit in Eastern Qinling Mountains (Central China) and its implication for estimating age of the pleistocene lithic artifacts. Quaternary Sci. 2007;27:559–67.

    Google Scholar 

  50. Sun SS, Guo YL, Favre A, Fu PC. Genetic differentiation and evolutionary history of two medicinal gentians (Gentiana Stipitata Edgew. And Gentiana Szechenyii Kanitz) in the Qinghai-Tibet Plateau. J Appl Res Med Aromat Plants. 2022;30:100375. https://doi.org/10.1016/j.jarmap.2022.100375.

    Article  CAS  Google Scholar 

  51. Fu PC, Sun SS, Khan G, Dong XX, Tan JZ, Favre A, et al. Population subdivision and hybridization in a species complex of Gentiana in the Qinghai-Tibetan Plateau. Ann Bot. 2020;125(4):677–90. https://doi.org/10.1093/aob/mcaa003.

    Article  PubMed  PubMed Central  Google Scholar 

  52. Fu PC, Favre A, Wang R, Huang Y, Sun SS. Between allopatry and secondary contact: differentiation and hybridization among three sympatric Gentiana species in the Qinghai-Tibet Plateau. BMC Plant Biol. 2022;22(1):504. https://doi.org/10.1186/s12870-022-03879-0.

    Article  PubMed  PubMed Central  Google Scholar 

  53. Fu PC, Twyford AD, Hao YT, Zhang Y, Chen SL, Sun SS. Hybridization and divergent climatic preferences drive divergence of two allopatric Gentiana species on the Qinghai–Tibet Plateau. Ann Bot. 2023;132:1271–88. https://doi.org/10.1093/aob/mcad179.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Wu SD, Wang Y, Wang ZF, Shrestha N, Liu JQ. Species divergence with gene flow and hybrid speciation on the Qinghai–Tibet Plateau. New Phytol. 2022;234:392–404. https://doi.org/10.1111/nph.17956.

    Article  CAS  PubMed  Google Scholar 

  55. Dai L, Guo X, Zhang F, Du Y, Ke X, Li Y, et al. Seasonal dynamics and controls of deep soil water infiltration in the seasonally-frozen region of the Qinghai-Tibet plateau. J Hydrol. 2019;571:740–8. https://doi.org/10.1016/j.jhydrol.2019.02.021.

    Article  Google Scholar 

  56. Zhang F, Li H, Li Y, Guo X, Dai L, Lin L, et al. Strong seasonal connectivity between shallow groundwater and soil frost in a humid alpine meadow, northeastern Qinghai-Tibetan Plateau. J Hydrol. 2019;574:926–35. https://doi.org/10.1016/j.jhydrol.2019.05.008.

    Article  Google Scholar 

  57. Rybczyński JJ, Davey MR, Mikuła A. The Gentianaceae-volume 2: biotechnology and applications. Heidelberg: Springer; 2015.

    Book  Google Scholar 

  58. Fu PC, Sun SS, Hollingsworth PM, Chen SL, Favre A, Twyford AD. Population genomics reveal deep divergence and strong geographical structure in gentians in the Hengduan Mountains. Front Plant Sci. 2022a;13:936761. https://doi.org/10.3389/fpls.2022.936761.

    Article  PubMed  PubMed Central  Google Scholar 

  59. Elshire RJ, Glaubitz JC, Sun Q, Albers CA, Banks E, DePristo MA, et al. A robust, simple genotyping-by-sequencing (GBS) approach for high diversity species. PLoS ONE. 2011;6:e19379. https://doi.org/10.1371/journal.pone.0019379.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  60. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114–20. https://doi.org/10.1093/bioinformatics/btu170.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  61. Li H, Durbin R. Fast and accurate short read alignment with Burrows–Wheeler transform. Bioinformatics. 2009;25:1754–60. https://doi.org/10.1093/bioinformatics/btp324.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  62. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25:2078–9. https://doi.org/10.1093/bioinformatics/btp352.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  63. Li T, Yu X, Ren YM, Kang MH, Yang WJ, Feng LD, et al. The chromosome-level genome assembly of Gentiana Dahurica (Gentianaceae) provides insights into gentiopicroside biosynthesis. DNA Res. 2022;29:dsac008. https://doi.org/10.1093/dnares/dsac008.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  64. Zhu MJ, Wang ZY, Yang YZ, Wang ZF, Mu WJ, Liu JQ. Multi-omics reveal differentiation and maintenance of dimorphic flowers in an alpine plant on the Qinghai-Tibet Plateau. Mol Ecol. 2023;32:1411–24. https://doi.org/10.1111/mec.16449.

    Article  CAS  PubMed  Google Scholar 

  65. Tarasov A, Vilella AJ, Cuppen E, Nijman IJ, Prins P. Sambamba: fast processing of NGS alignment formats. Bioinformatics. 2015;31:2032–4. https://doi.org/10.1093/bioinformatics/btv098.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  66. Garrison E, Marth G. Haplotype-based variant detection from short-read sequencing. arXiv Preprint arXiv. 2012. https://doi.org/10.48550/arXiv.1207.3907. 1207.3907.

    Article  Google Scholar 

  67. Danecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA, et al. The variant call format and VCFtools. Bioinformatics. 2011;27:2156–8. https://doi.org/10.1093/bioinformatics/btr330.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  68. Goudet J. Hierfstat, a package for r to compute and test hierarchical F-statistics. Mol Ecol Notes. 2005;5:184–6. https://doi.org/10.1111/j.1471-8286.2004.00828.x.

    Article  Google Scholar 

  69. Weir BS, Cockerham CC. Estimating F-statistics for the analysis of population structure. Evolution. 1984;38:1358–70. https://doi.org/10.2307/2408641.

    Article  CAS  PubMed  Google Scholar 

  70. Raj A, Stephens M, Pritchard JK. fastSTRUCTURE: variational inference of population structure in large SNP data sets. Genetics. 2014;197:573–89.

    Article  PubMed  PubMed Central  Google Scholar 

  71. Rosenberg NE. DISTRUCT: a program for the graphical display of population structure. Mol Ecol Notes. 2004;4:137–8. https://doi.org/10.1046/j.1471-8286.2003.00566.x.

    Article  Google Scholar 

  72. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MA, Bender D, et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. Amer J Hum Genet. 2007;81(3):559–75. https://doi.org/10.1086/519795.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  73. Hoang DT, Chernomor O, Von Haeseler A, Minh BQ, Vinh LS. UFBoot2: improving the ultrafast bootstrap approximation. Mol Biol Evol. 2018;35:518–22. https://doi.org/10.1093/molbev/msx281.

    Article  CAS  PubMed  Google Scholar 

  74. Nguyen LT, Schmidt HA, Von Haeseler A, Minh BQ. IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. Mol Biol Evol. 2015;32:268–74. https://doi.org/10.1093/molbev/msu300.

    Article  CAS  PubMed  Google Scholar 

  75. Favre A, Pringle JS, Heckenhauer J, Kozuharova E, Gao QB, Lemmon EM, et al. Phylogenetic relationships and sectional delineation within Gentiana (Gentianaceae). Taxon. 2020;69:1221–38. https://doi.org/10.1002/tax.12405.

    Article  Google Scholar 

  76. Chen CL, Zhang L, Li JL, Mao XX, Zhang LS, Hu QJ, et al. Phylotranscriptomics reveals extensive gene duplication in the subtribe Gentianinae (Gentianaceae). J Syst Evol. 2021;59:1198–208. https://doi.org/10.1111/jse.12651.

    Article  Google Scholar 

  77. Fu PC, Sun SS, Twyford AD, Li ZH, Wang XM, et al. Lineage-specific plastid degradation in subtribe Gentianinae (Gentianaceae). Ecol. Evol. 2021;11:3286–99. https://doi.org/10.1002/ece3.7281.

    Article  Google Scholar 

  78. Drummond AJ, Suchard MA, Xie D, Rambaut A. Bayesian phylogenetics with BEAUti and the BEAST 1.7. Mol Biol Evol. 2012;29:1969–73. https://doi.org/10.1093/molbev/mss075.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  79. Bouckaert R, Heled J, Kühnert D, Vaughan T, Wu CH, Xie D, et al. BEAST 2: a software platform for bayesian evolutionary analysis. PLoS Comput Biol. 2014;10:e1003537. https://doi.org/10.1371/journal.pcbi.1003537.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  80. Mai DH. Die mittelmiozänen und obermiozänen Floren Aus Der Meuroer Und Raunoer Folge in Der Lausitz. Teil I: Farnpflanzen, Koniferen Und Monokotyledonen. Palaeontographica Abteilung B Pälaophytologie. 2000;256:1–68.

    Article  Google Scholar 

  81. Schenk JJ, Axel J. Consequences of secondary calibrations on divergence time estimates. PLoS ONE. 2016;11:e0148228. https://doi.org/10.1371/journal.pone.0148228.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  82. Janssens S, Couvreur TLP, Mertens A, Dauby G, Dagallier LP, Vanden Abeele SV, et al. A large-scale species level dated angiosperm phylogeny for evolutionary and ecological analyses. Biodivers Data J. 2020;8:e39677. https://doi.org/10.3897/BDJ.8.e39677.

    Article  PubMed  PubMed Central  Google Scholar 

  83. Collin F, Durif G, Raynal L, Lombaert E, Gautier M, Vitalis R, et al. Extending approximate bayesian computation with supervised machine learning to infer demographic history from genetic polymorphisms using DIYABC random forest. Mol Ecol Resour. 2021;21:2598–613. https://doi.org/10.1111/1755-0998.13413.

    Article  PubMed  PubMed Central  Google Scholar 

  84. Fitzpatrick MC, Chhatre VE, Soolanayakanahally RY, Funch P. Experimental support for genomic prediction of climate maladaptation using the machine learning approach gradient forests. Mol Ecol Resour. 2021;21:2749–65. https://doi.org/10.1111/1755-0998.13374.

    Article  CAS  PubMed  Google Scholar 

  85. Ellis N, Smith SJ, Pitcher CR. Gradient forests: calculating importance gradients on physical predictors. Ecology. 2012;93:156–68. https://doi.org/10.1890/11-0252.1.

    Article  PubMed  Google Scholar 

  86. Dixo P. VEGAN, a package of R functions for community ecology. J Veg Sci. 2003;14:927–30. https://doi.org/10.1111/j.1654-1103.2003.tb02228.x.

    Article  Google Scholar 

  87. Hijmans RJ, Van Etten J, Sumner M, Cheng J, Baston D, Bevan A et al. Package ‘raster’. R package version. 2015; 3.6–26. https://cran.r-project.org/web/packages/raster/index.html (8 October 2023, date last accessed).

  88. Goslee SC, Urban DL. The ecodist package for dissimilarity-based analysis of ecological data. J Stat Softw. 2007;22:1–19. https://doi.org/10.18637/jss.v022.i07.

    Article  Google Scholar 

  89. Mantel N. The detection of disease clustering and a generalized regression approach. Cancer Res. 1967;27:209–20.

    CAS  PubMed  Google Scholar 

Download references

Acknowledgements

We thank Fei Liu to offer help in statistics analyses, and thank Ming Tang for sample collection from his project “The Survey of Wildlife Resources in Key Areas of Xizang (Phase II) (ZL202303601)”.

Author information

Authors and Affiliations

Authors

Contributions

P.C.F. designed the experiments, collected the samples, analyzed the data and wrote the main manuscript; B.J.M, H.X.W. and S.W.Y. analyzed the data; R.X. collected the samples; S.S.S. designed the experiments, prepared the figures and revised the manuscript.

Corresponding author

Correspondence to Shan-Shan Sun.

Ethics declarations

Ethics approval and consent to participate

We comply with the IUCN Policy Statement on Research Involving Species at Risk of Extinction and the Convention on the Trade in Endangered Species of Wild Fauna and Flora and confirm that all methods were performed in accordance with the relevant guidelines/regulations/legislation. These wild species samples in this study have not been included in the list of key protected plants. The sampling of wild species in this study was met local policy requirements and did not affect the survival or reproduction of the species. The materials were identified by Dr. Pengcheng Fu, and their voucher specimens were deposited in the herbarium of Luoyang Normal University.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Electronic supplementary material

Below is the link to the electronic supplementary material.

Supplementary Material 1

Supplementary Material 2

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Fu, PC., Mo, BJ., Wan, HX. et al. Divergence of alpine plant populations of three Gentianaceae species in the Qinling sky Island. BMC Plant Biol 25, 144 (2025). https://doi.org/10.1186/s12870-025-06165-x

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12870-025-06165-x

Keywords