- Research article
- Open Access
Nunataks or massif de refuge? A phylogeographic study of Rhodiola crenulata (Crassulaceae) on the world’s highest sky islands
BMC Evolutionary Biology volume 18, Article number: 154 (2018)
Quaternary climatic oscillations had tremendous effects on the current distribution of species. Here, we aim to elucidate the glacial history of Rhodiola crenulata, a perennial herb almost exclusively restricted to rock crevices on mountain peaks, and to test whether the nunatak or massif de refuge hypotheses could explain its distribution pattern.
Six haplotypes and six ribotypes were detected in the cpDNA data set and the ITS data set, respectively. The divergence of R. crenulata and its closest relatives was dated have occurred ca. 0.65 Mya, during the Naynayxungla glaciation on the QTP. Mismatch distribution analysis suggested that the species experienced a range expansion around 0.31 Mya. Populations with high genetic and haplotype diversity were found on the QTP platform as well in the Hengduan Mountains. The ecological niche modeling results showed that there were suitable habitats on both the QTP platform and in the Hengduan Mountains during the LGM.
Our results support a scenario that both nunataks and the massif de refuge hypotheses could explain the distribution of R. crenulata. We also confirmed that Quaternary climatic oscillations could promote plant speciation in some circumstances. This study adds to a growing body of evidence suggesting that the QTP plant lineages exhibited diverse reactions to the Quaternary climatic oscillations.
Mountains have long been recognized as island-like systems [1,2,3]. These so-called “sky islands” could provide us with important insights into evolutionary mechanisms of island systems. Limited size, distinct boundaries, isolation and dispersal limitation are shared attributes of these systems. Researchers since Darwin have documented the divergence of populations in relation to allopatry, evolutionary modification of form and function, and strikingly rapid diversification on island systems, including sky islands both on the Qinghai-Tibetan Plateau (QTP) and in the Andes . However, how plants on such sky islands responded to the Quaternary climate oscillations, which had tremendous effects on species ranges, has not been fully explored.
The current distribution of genetic lineages of species in the Northern Hemisphere was greatly influenced by the Quaternary climatic oscillation [5, 6]. During the ice ages, the decreasing temperature and reduced altitude of the mountain snow line  allowed organisms originally living in alpine zones either to move to lower elevations or to find a suitable refugium in situ. Two hypotheses exist relating to how alpine organisms survive glacial periods based on studies in the Alps and the north Atlantic area [8,9,10]: the “nunatak” hypothesis suggests survival on mountain tops within the ice sheet , whilst the massif de refuge hypothesis suggests that the periphery of mountain ranges provided large refugia for alpine species, allowing the refugees to re-occupy alpine zones after the retreat of glaciers . Previous phylogeographic studies conducted on the QTP, one of the largest and highest plateaus in the world, provide us with important insights into how plants on the plateau reacted to the Pleistocene climate oscillation. Three patterns have been revealed: contraction to the Hengduan Mountains during glaciations and subsequent recolonization onto the plateau [13, 14], which corresponds to the massif de refuge hypothesis; surviving in refugia that existed on the plateau platform during glaciations and local expansion during interglacial or postglacial [15, 16]; and the “microrefugia” hypothesis (i.e., multiple microrefugia during the Last Glacial Maximum (LGM, ca. 20 kya) across the whole of the current distribution range of species), which was confirmed by two previous studies of woody plants on the QTP [17, 18]. Although such studies shed important light on the evolutionary history of alpine species on the QTP during the glacial period, most of them have focused on woody plants or herbs of alpine meadows. Few studies have been conducted to illustrate how species growing in rock crevices on mountain peaks reacted to the Quaternary climatic oscillations.
Rhodiola crenulata (HK. f. et Thoms.) H. Ohba, our study species, is almost exclusively restricted to rock crevices on mountain peaks (Fig. 1), from 4300 m to 5600 m in elevation , which makes it one of the highest living vascular plants on the QTP. Such a distribution pattern represents a classical continental island system, with mountain peaks regarded as isolated islands. The species’ affinity for high alpine and nival altitudes suggests that it would survive on nunataks, which certainly offered extremely harsh conditions for survival. Consequently, our first goal was to elucidate the glacial history of R. crenulata, and to test the nunatak and massif de refuge hypotheses. Rhodiola crenulata has been used as an important source of adaptogens, hemostatics, and tonics in traditional Tibetan medicines for thousands of years . The accelerated and uncontrolled exploitation of R. crenulata has made it extremely endangered. Thus our second goal was to explore the genetic diversity and distribution pattern of this species to offer more information for conservation and management programs. To achieve our goals, we examined plastid and nuclear ITS markers and used ecological niche models (ENMs) to conduct a phylogeographic study. The ENMs that are widely used in phylogeographic studies, can provide important evolutionary and biological insights by allowing the evaluation of biogeographic hypotheses or to interpret genetic diversity patterns [21, 22].
Our sampling did not require specific permits for the locations involved as we held a permission letter from College of Life Sciences, Shaanxi Normal University, Xi’an relating to colleting samples. We sampled 16 populations across the distribution range of R. crenulata during the summers of 2016 to 2017. For each population, we generally sampled 10–20 individuals. Each sampled individual in every population was at least 20 m away from any other to avoid sampling clonal individuals. In total, we collected 253 individuals and dried their leaves with silica gel. Based on a previous phylogenetic study, we chose one individual of R. wallichiana (Hk.) S. H. Fu to act as an outgroup .
DNA extraction, PCR amplification, cloning and sequencing
We used a Plant Genomic DNA Kit (TianGen Biotech, Beijing, China) to extract DNA from silica-gel dried leaves. PCR amplification was conducted in a 20 μl system (2 μl 10 × buffer, 0.5 μl of each primer, 0.4 μl of dNTP mixture, 1 U of Taq polymerase, 1 μl template genomic DNA), using primers for ITS, ITS-1 and ITS-4 , for trnL-F, c and f  and for trnS-G, trnS and trnG . We used the PCR conditions described by Liu et al. . For chloroplast sequences and most of the ITS sequences, direct sequencing was conducted. For those ITS sequences with multiple peaks, we ligated the PCR product into pGEM-T Easy Vector using a Promega Kit (Promega Corporation, Madison, WI, USA), and chose a plasmid containing the right insertion to sequence. We conducted all sequencing on a 3730 automatic DNA sequencer at Tsingke Biotech, Xi’an, China. The contigs obtained were edited and assembled using ContigExpress (a component of Vector NTI Suite 6.0, InforMax). ClustalW v. 1.7  was used to align the data set, and BioEdit v. 7.0.1. was used to check the aligned sequence by eye.
Phylogenetic analysis and divergence time inference
Sequences from two cpDNA regions (trnL-F and trnS-G) were concatenated into a matrix because the chloroplast genome can be considered to have evolved as a single entity in plants. Segregating sites and haplotypes were extracted using DnaSP v5 . To infer phylogenetic relationships between extracted haplotypes/ribotypes, we applied the Maximum Parsimony method in PAUP* v. 4.0b10  and then conducted Bayesian inference analysis. We specified the following parameters for the parsimony calculation: heuristic searches with MULTREES, TBR branch swapping, 1000 replicates of random addition sequences for starting trees. Bootstrap values  were obtained from a run of 1000 replicates, with ten replicates of random addition sequences and NNI branch swapping.
Because Bayesian inference requires a nucleotide substitution model, we determined one for each data set with the Akaike information criterion (AIC) in Modeltest v. 3.7 [32, 33] (TPM2uf + G for the plastid data set, and SYM + G for the ITS data set). MrBayes v. 3.2.1 (nst = 6, rates = gamma for the plastid data set, nst = 6, rates = gamma for the ITS data set)  was used the inference. We surveyed 10,000,000 generations of the four Metropolis-coupled Markov chains, and sampled every 1000 generations. The average stand deviation of split frequencies was used to assess the convergence of two independent runs, and the threshold was set to 0.05. We discarded the first 25% of generations as burn-in, and constructed a 50%-majority rule consensus tree with the remainder. The posterior proportion for each node was used to estimate robustness of the BI trees. In addition, we constructed haplotype/ribotype networks with NETWORK v. 126.96.36.199  to reveal relationships between sequences with shallow genetic divergences.
To infer the divergence time between lineages, we first examined the molecular clock hypothesis with a likelihood ratio test in PAUP* v. 4.0b10. The results showed that both the ITS and cpDNA data sets fit the molecular clock hypothesis (2logeLR = 1.36, df = 5, p > 0.05 for the ITS data, and 2logeLR = 0.78, df = 5, p > 0.05 for the plastid data). We then used both the ITS and plastid data sets to conduct a dating analysis using the BEAST software . All parameters were sampled every 1000 generations from 10,000,000 MCMC generations with the first 25% as burn-in. Tracer  was used to examine the convergence of chains. As far as we know, there are no reliable fossils of R. crenulata and its relatives available; we therefore used substitution rates to estimate divergence times. It has been reported that the substitution rate of nrITS in shrubs and herbs varies from 3.46 × 10− 9 to 8.69 × 10− 9 s/s/y , thus we used a normal distribution prior to cover these ranges with a 95% confidence interval (6.075 × 10− 9 1.590 × 10− 9). The evolutionary rate of the plastid markers was assumed to be 2 × 10− 9 s/s/y .
We used DnaSP v5  to estimate the population genetic diversity indexes: haplotype diversity (h) and nucleotide diversity (π) for each population (hS, πS) and at the species level (hT, πT). We also calculated average gene diversity within the population (HS), total gene diversity (HT) and between population divergence (GST, NST) using the program PERMUT (available at https://www6.bordeaux-aquitaine.inra.fr/biogeco/Production-scientifique/Logiciels/Contrib-Permut/Permut) with 1000 permutation tests. If the NST value is larger than GST, it means that closely related haplotypes will appear in the same area with higher probability than remotely related ones, thus indicating the presence of phylogeographic structure . The Mantel test with 1000 random permutations was implemented to test the significance of isolation by distance between populations using ARLEQUIN v. 3.5 .
To see if there is any spatial structure in populations of R. crenulata, we used SAMOVA v. 1.0  to assess the spatial genetic pattern. By maximizing the FCT value, the program can search different K values that the user defines as groups of geographically adjacent populations. We ran the program repeatedly with K values from 2 to 10. After defining geographic groups, we calculated the amount of variation among populations within a geographic group and within a population using the hierarchical analysis of molecular variance (AMOVA) framework, which was run using ARLEQUIN v. 3.5 . We used a nonparametric permutation procedure with 1000 permutations to test for any significant difference.
Then we inferred signatures of demographic expansion by calculating Tajima’s D and Fu’s Fs statistics [42, 43] and conducting a mismatch distribution analysis [44, 45]. D and Fs statistics are significantly negative if there is population expansion, as an excess of rare new mutations will originate during the expansion process. A previous study has shown that mismatch distribution is not affected by population structure , thus we pooled haplotypes of each clade together to conduct the analysis. The difference between observed mismatch distributions and expected distribution under a sudden expansion model [44, 47] was tested by examining the sum of squared deviations (SSD) and the raggedness index (HRag)  with 1000 parametric bootstrap replicates. If expansion was detected in a group, the parameter-value τ was used to estimate the expansion time (in generations) with the equation t = τ/2u [44, 46], where u = μ × k × g, in which μ represents substitution rate (s/s/y), k represents the length of aligned matrix, and g is the generation time in years. For the present study, k was 1531 bp, and the substitution rate for the cpDNA genome was assumed to be 2 × 10− 9 s/s/y [38, 49]. Generation time was estimated to be 10 years based on personal observation.
Ecological niche modeling
We used ENMs to examine the potential range shifts of R. crenulata in response to glacial climatic oscillations. Locations of R. crunulata were obtained from field collections of the authors and from herbarium records online (e.g., Chinese Virtual Herbarium, National Specimen Information Infrastructure, and Global Biodiversity Information Facility). Data from online databases were checked to exclude misidentification. A total of 88 spatially unique localities were used for analysis. MAXENT v. 3.3.3e  was used to model the ecological niches of R. crenulata. Twenty replicates were employed, from which 80% of the distribution coordinates were used for training and 20% for testing. Environmental layers of 19 bioclimatic variables (Additional file 1: Table S1) for the Last Glacial Maximum (LGM), Last Interglacial (LIG), and the current time and the near future (year 2080) were downloaded from the WorldClim dataset at a spatial resolution of 2.5 arc-minutes and these were employed for the modeling . To exclude highly correlated climate variables, pairwise correlations were examined among the 19 variables within the distribution area of R. crenulata. The nine variables (Additional file 1: Table S1) with pairwise Pearson correlation coefficients below 0.7 were used. Area under the “receiver operating characteristic (ROC) curve” (AUC) [52, 53] values were used to evaluate the accuracy of each model prediction. The threshold for good performance was set to 0.7 . DIVA-GIS v. 7.5  was used to draw the range of suitable distributions.
Plastid DNA sequence variation and distribution
The total alignment length of the two plastid DNA sequences was 1531 bp. In combination, we revealed six haplotypes determined by seven nucleotide substitutions. Of the 16 populations sampled, 13 were fixed for a single haplotype, and the remaining three were polymorphic (Table 1 & Additional file 2: Table S2; Fig. 2). Population ML harbored three haplotypes, which was the most of any of the populations. Three haplotypes (H4, H5 and H6) occurred only in one population each, while H2 was shared by 12 of the 16 populations. At species level, haplotype diversity was h = 0.554. Considering the populations, the highest haplotype diversity was detected in the ML population, with a value of 0.353 (Table 1). Nucleotide diversity was pi = 0.00063 at the overall scale, and the highest value was found in population QE-3, with 0.00068 (Table 1). Total gene diversity (HT) was much higher than within-population gene diversity (HS) (0.598 and 0.049, respectively, Table 2).
NrDNA ITS variation and ribotype distribution
The total length of aligned ITS sequences was 606 bp, with a length variation from 604 to 606 bp. We detected six different ITS sequences (ribotypes) in all 253 individuals, which were determined by four substitutions. Among the six ribotypes, three occurred in only one population. R1 was the most widespread ribotype, occurring in 15 of the 16 populations we sampled (Table 1 & Additional file 3: Table S3; Fig. 3). Ribotype diversity was h = 0.401 at the species level. As with the cpDNA data, population ML had the highest h values (Table 1). The total nucleotide diversity was estimated to be 0.0007, and considering the populations, the highest was detected in the ML population (Table 1). Within-population gene diversity (HS) was also much lower than total gene diversity (HT) (0.142 and 0.375, respectively; Table 2) as demonstrated in the cpDNA data. Strikingly, 11 out of the 16 sampled populations only harbored a single ribotype.
Phylogenetic relationships and lineage divergence
For both data sets, MP and Bayesian inference yielded largely congruent tree topologies, thus we only show the Bayesian tree with the MP bootstrap value marked on each branch. The phylogenetic relationship between haplotypes is shown in Fig. 4a. H3 and H6 comprised a moderately-supported clade (Fig. 4a), and together with four other haplotypes formed a polytomy. The haplotype network (Fig. 2b) revealed a consistent relationship between detected haplotypes within the phylogenetic tree (Fig. 4a),
Phylogenetic relationships for ITS ribotypes are shown in Fig. 4b. Ribotype 2 and other ribotypes diverged first, while the remaining five ribotypes formed a clade (BP = 63, PP = 0.99). The network relationships of ribotypes are shown in Fig. 3b, which is consistent with the phylogenetic trees (Fig. 3a), but exhibits a reticulate relationship.
The dating analysis based on the ITS data showed that R. crenulata diverged from its closest relatives 0.65 Mya (95% HPD: 0.20–1.20 Mya), i.e. in the Pleistocene (Fig. 5). Further divergence of R2 and the other five ribotypes occurred 0.41 Mya (95% HPD: 0.13–0.76 Mya). Dating based on the cpDNA data set was consistent with that based on ITS data but with earlier dates: the divergence of R. crenulata and R. wallichiana occurred 1.25 Mya (95% HPD: 0.44–2.22 Mya), still in the Pleistocene (Additional file 4: Figure S1) and further divergence between H3, H6 and other haplotypes occurred 0.82 Mya (95% HPD: 0.29–1.44 Mya).
The GST values (0.918 and 0.622 for cpDNA and ITS data, respectively) (Table 2) showed high differentiation between R. crenulata populations. However, the permutation tests for both cpDNA and ITS data sets showed that NST was not significantly higher than GST (Table 2). In the SAMOVA analysis for the cpDNA data set, the FCT value reached a plateau with K = 4 (Additional file 5: Figure S2). Population ML formed a group on its own, and DD-1 and DD-2 formed another group. The remaining populations were separated into two groups: JCL, SJL-1 and SJL-2 were in one group and the others in another group. For the ITS data set, K = 2 resulted in the maximum value, thus two groups were detected. JCL was an independent group, and the remaining populations formed another group. AMOVA analysis showed that for the cpDNA data set, 93.91% of the total variations occurred among groups according to SAMOVA. For the ITS data set, variations among groups also accounted for 69.72% of the total variation (Table 3). We found no pattern of geographic isolation revealed by the Mantel test based on both data sets.
Under the spatial population expansion model, a unimodel mismatch distribution was detected including all haplotypes (Additional file 6: Figure S3). We also found that observed variance (SSD) and the raggedness index were not significantly different from the expected values (Table 4). In addition, a negative value for Tajima’s D was calculated (Table 4). Thus, all evidence suggested that R. crenulata underwent a rapid expansion, the time of which was dated to around 0.31 Mya (Table 4). However, we failed to detect population expansion signals in the individual clades shown in Fig. 4a.
Species distribution models for R. crenulata
ENMs performed with high predictive ability, with AUC > 0.98 and standard deviation [SD] < 0.01 in each model. The predicted potential distributions of R. crenulata during the LIG, LGM, present day, and in the future are shown in Fig. 6. Our current sampling area covers most of the current potential distribution area of R. crenulata (Fig. 1c). In both the LIG and LGM models, the distribution area was smaller and also more fragmented than at present. Interestingly, the LIG model showed a more significant contraction to the Hengduan Mountain area (Fig. 6a). Focusing on higher habitat suitability (> 0.5), the existence of suitable habitats on the QTP was predicted in the LGM model (Fig. 6b), but not in the LIG model. We also modeled the potential distribution of R. crenulata under future climate change scenarios (year 2080) (Fig. 6d). We found that in 2080, R. crenulata will expand into the Himalaya area slightly, but will contract its distribution in the Hengduan Mountains compared to the present.
Phylogeographic pattern of R. crenulata
Many authors have investigated how the Quaternary climate oscillations have shaped the distributions patterns of plants on the QTP and its adjacent areas [15,16,17,18, 39, 56,57,58]. Some species appear to comply with the massif de refuge hypothesis : Juniperus przewalskii  and Picea crassifolia  showed reduced genetic diversity or even single haplotypes on the QTP platform but high levels of genetic diversity and private haplotypes in the Hengduan Mountains. On the other hand, refugia existed both on the plateau and in the Hengduan Mountains for Potentilla glabra  and Aconitum gymnandrum . In addition, the Juniperus tibetica complex  and Hippophae tibetana  provide evidence that “microrefugia” existed across the distributional range.
Our results represent an interesting case: neither the nunataks nor the massif de refuge hypotheses can be rejected for R. crenulata, which is also the case for Potentilla glabra  and Aconitum gymnandrum . Glacial refugia are characterized by high levels of genetic diversity and the presence of private haplotypes . In the case of R. crenulata, both cpDNA and ITS data sets showed high genetic diversity and private haplotypes in population ML (Figs. 2 and 3; Table 1), which is on the plateau, and in QE and HS populations (Figs. 2 and 3; Table 1), which are located in the Hengduan Mountains. The ecological modeling results also showed that, during the LGM, there were suitable habitats for R. crenulata both on the QTP and in the Hengduan Mountains (Fig. 6b). These results indicate that glacial refugia exited both on the QTP platform and in the Hengduan Mountains.
In previous studies, it has been suggested that the south-eastern edge of the QTP and the Hengduan Mountians represented important refugia for the QTP species [16, 23, 56]. In the cpDNA data, populations located in this area, such as DD-1 and QE-3 had high genetic diversity and haplotype diversity (Table 1). In addition, population QE-3 harbored a private haplotype H6 (Fig. 2). A similar pattern can be seen in the ITS data set: population QE-1, XC and HS-1 had high genetic diversity and high haplotype diversity. Population QE-1 harbored private ribotype R5. Thus, our results confirmed that the south-eastern edge of the QTP and the Hengduan Mountains played an important role in hosting the genetic diversity of QTP plant species. Another population with high genetic diversity and haplotype uniqueness is ML, which is on the QTP. The refugium status was supported by both cpDNA and the ITS datasets. Interestingly, a previous study that focused on another Rhodiola species, R. alsia, also revealed that this area was a potential refugium for the species. Similarly, studies on Potentilla glabra  and Aconitum gymnandrum  also revealed in situ refugia on the QTP. Rhodiola crenulata often grows on schist on mountain slopes, rocky places and rock crevices at 4300–5600 m in elevation . It is reasonable to suppose that it would have survived in situ on the QTP on nunataks that were ice free.
Our AMOVA analysis results suggest that 93.9% and 69.7% of the genetic diversity was found among groups defined by the SAMOVA analysis for the cpDNA and ITS data, respectively. It is generally expected that such genetic differentiation should be coupled with a distinct phylogeographic structure . However, we found that neither cpDNA nor ITS datasets showed significantly higher NST than GST (Table 2), which indicated no phylogeographic structure across the distribution area of R. crenulata. We note that three of the six haplotypes (50%) and four out of the six ribotypes (67%) are endemic to one population. Furthermore, only one haplotype and one ribotype were widespread across the distribution range (Figs. 2 and 3). A small proportion of the widespread haplotype or ribotype and a larger proportion of endemic ones may cause the absence of phylogeographic structure in R. crenulata. Notably, H2 appeared in 12 of the 16 populations sampled, and R1 in 15 of the 16 populations. This pattern may have been the result of recent range expansion from glacial refugia. Mismatch distribution analysis based on the cpDNA dataset showed that all populations of R. crenulata as a whole experienced a range expansion, as indicated by negative Tajima’s D and a unimodal mismatch distribution (Table 4; Additional file 5: Figure S2). This expansion event was dated to 0.31 Mya, at the time of the interglacial between the Naynayxungla Glaciation (0.5–0.72 Mya) and the Guxiang Glaciation (0.13–0.3 Mya) . The wide distribution of H2 and R1 may have been caused by this expansion event.
We also found a very low genetic diversity compared to congeneric species [23, 57] and other species in the same region [13,14,15,16, 39]. According to our results, this extremely low genetic diversity may be due to the species’ recent origin, as well as low within population diversity. In our field study, populations of R. crenulata are narrowly confined to mountain peaks, with a limited population size. This may be the reason for the low within population genetic diversity. Another possibility is that “Orbitally Forced Range Dynamics” – ORD [62, 63] have resulted in the current pattern of limited haplotype diversity. As R. crenulata is confined to glacial relic rocks, the Milankovitch climate oscillations may have repeatedly wiped out new mutations, as individuals carrying these mutations often went extinct during these oscillations.
Implications for quaternary speciation
Although the role that the Pleistocene played in shaping intraspecific genetic diversity has been examined by a series of studies [5, 6, 64, 65], it is still debated whether speciation could have been triggered by the Pleistocene climate oscillations . Palaeontological evidence has shown that most species remained unchanged throughout the Pleistocene . However, an increasing number of studies reported speciation in the Quaternary [67,68,69, 70] and, indeed, the role of Quaternary speciation has attracted attention since Darwin . Our dating analysis based on both cpDNA and ITS data sets clearly demonstrated that divergence between R. crenulata and its closest relative R. wallichiana occurred about 0.65 Mya, during the Naynayxungla glaciation on the QTP . Although absolute dating of molecular phylogenies using the molecular clock method needs to be treated cautiously, our results suggest that the speciation of R. crenulata took place during the Pleistocene. The cpDNA data also suggested the divergence between R. crenulata and R. wallichiana be in the Pleistocene (Additional file 4: Figure S1). The evidence is in accordance with the fact that R. crenulata only grows on glacial relics, i.e. schist on mountain slopes, rocky places and rock crevices at 4300–5600 m in elevation. There were no glacial relics before the Pleistocene, so there were no suitable habitats for R. crenulata. It is reasonable to hypothesize that by adapting rapidly to glacial relics, R. crenulata diverged from its closest relative R. wallichiana, as the latter often grows at the moist margin of forests. Further studies based on extensive low copy nuclear markers and coalescence calculations are needed to elucidate the detailed evolutionary history of the two species.
Conservation and management in the future
The roots and rhizomes of R. crenulata have been included in the Pharmacopoeia of China  as the authentic Rhodiolae Crenulatae Radix et Rhizoma, where they are listed as enhancing inner spiritual power, concentration and physical endurance. Since the 1980s, the accelerated and uncontrolled use of R. crenulata in China has severely reduced its population. Population genetic studies with inter-simple sequence repeat markers have demonstrated low genetic diversity of this species . To preserve this threatened and endemic species, conservation programs should be planned and carried out. An ideal conservation program should involve all known populations, and populations that harbor higher genetic diversity, higher heterozygosity and unique haplotypes should be given priority. As our AMOVA analysis showed that most genetic variation exists among groups, conservation programs should be undertaken in each of the four groups. In particular, populations on Mila Mt. and Queer Mt. should be given the highest priority, because populations there had the highest genetic diversity and harbored multiple unique haplotypes (Figs. 2 and 3).
Another suggestion came from our ENMs analysis, which is helpful when developing conservation and management strategies because it predicts potential changes in geographic distribution under future climate change [73, 74]. It is reasonable to hypothesize that in a warming climate, mountain species will migrate toward higher altitudes . As R. crenulata already grows at very high altitudes, as the climate becomes warmer, its suitable habitats will shrink. This hypothesis was confirmed by our ENMs (Fig. 6d): compared to the current distribution, the distribution of R. crenulata would shrink into two major areas, one in the Hengduan Mountains and the other one on the QTP. In addition, our ENMs suggested that based on the global warming model (Fig. 6d), populations from lower altitudes will experience lower habitat suitability (< 0.5). From the perspective of conservation management, populations that face a higher extinction risk, which are in low altitudes, should be given priority for protection (e.g., XL, DQ). Areas around Mila Mt. and Queer Mt. all possess high habitat suitability (> 0.5), therefore should be listed as the core area for conservation of this species; this is supported by both our genetic data and ENMs.
We used two molecular markers and ecological niche modeling to elucidate the glacial history of R. crenulata, a perennial herb almost exclusively restricted to rock crevices on mountain peaks. Our results detected a strikingly recent divergence of R. crenulata from its closest relative, which was dated to around 0.65 Mya, during the Naynayxungla glaciation on the QTP. This result suggested that Quaternary climatic oscillations may have promoted plant speciation. In addition, we found that refugia existed in both the Hengduan Mountain area and on the QTP platform during the LGM, demonstrating a scenario that both the nunataks and the massif de refuge hypotheses apply to R. crenulata. Our study adds to a growing body of evidence suggesting that the QTP plant lineages exhibited diverse reactions to the Quaternary climatic oscillations.
Akaike Information Criterion
Analysis of Molecular Variance
Area Under the Curve
Ecological Niche Models
- HRag :
Raggedness index of Harpending
Last Glacial Maximum
Polymerase Chain Reaction
Receiver Operating Characteristic
Spatial Analysis of Molecular Variance
Sum of Squared Deviations
Hughes C, Eastwood R. Island radiation on a continental scale: exceptional rates of plant diversification after uplift of the Andes. Proc Natl Acad Sci U S A. 2006;103:10334–9.
Sklenář P, Hedberg I, Cleef AM. Island biogeography of tropical alpine floras. J Biogeogr. 2014;41:287–97.
Hughes CE, Atchison GW. The ubiquity of alpine plant radiations: from the Andes to the Hengduan Mountains. New Phytol. 2015;207:275–82.
Losos JB, Ricklefs RE. Adaptation and diversification on islands. Nature. 2009;457:830–6.
Hewitt G. The genetic legacy of the quaternary ice ages. Nature. 2000;405:907–13.
Hewitt GM. Genetic consequences of climatic oscillations in the quaternary. Philos Trans R Soc Lond Ser B Biol Sci. 2004;359:183–95.
Yugo O. Last glacial snowline altitude and paleoclimate of the eastern Asia. Quat Res. 1988;26:271–80.
Schoanswetter P, Tribsch A, Stehlik I, Niklfeld H. Glacial history of high alpine Ranunculus glacialis (Ranunculaceae) in the European Alps in a comparative phylogeographical context. Biol J Linn Soc. 2004;81:183–95.
Schneeweiss GM, Schonswetter P. A re-appraisal of nunatak survival in arctic-alpine phylogeography. Mol Ecol. 2011;20:190–2.
Wachter GA, Papadopoulou A, Muster C, Arthofer W, Knowles LL, Steiner FM, Schlick-Steiner BC. Glacial refugia, recolonization patterns and diversification forces in alpine-endemic Megabunus harvestmen. Mol Ecol. 2016;25:2904–19.
Dahl E. The nunatak theory reconsidered. Ecol Bull. 1987;38:77–94.
Nordal I. Tabula rasa after all? Botanical evidence for ice-free refugia in Scandinavia reviewed. J Biogeogr. 1987;14:377–88.
Yang FS, Li YF, Ding X, Wang XQ. Extensive population expansion of Pedicularis longiflora (Orobanchaceae) on the Qinghai-Tibetan plateau and its correlation with the quaternary climate change. Mol Ecol. 2008;17:5135–45.
Zhang YH, Volis S, Sun H. Chloroplast phylogeny and phylogeography of Stellera chamaejasme on the Qinghai-Tibet plateau and in adjacent regions. Mol Phylogenet Evol. 2010;57:1162–72.
Wang LY, Ikeda H, Liu TL, Wang YJ, Liu JQ. Repeated range expansion and glacial endurance of Potentilla glabra (Rosaceae) in the Qinghai-Tibetan plateau. J Integr Plant Biol. 2009;51:698–706.
Wang LY, Abbott RJ, Zheng W, Chen P, Wang YJ, Liu JQ. History and evolution of alpine plants endemic to the Qinghai-Tibetan plateau: Aconitum gymnandrum (Ranunculaceae). Mol Ecol. 2009;18:709–21.
Opgenoorth L, Vendramin GG, Mao K, Miehe G, Miehe S, Liu JL, Ziegenhagen B. Tree endurance on the Tibetan plateau marks the world's highest known tree line of the last glacial maximum. New Phytol. 2010;186:780.
Wang H, Qiong L, Sun K, Lu F, Wang Y, Song Z, Wu Q, Chen J, Zhang W. Phylogeographic structure of Hippophae tibetana (Elaeagnaceae) highlights the highest microrefugia and the rapid uplift of the Qinghai-Tibetan plateau. Mol Ecol. 2010;19:2964–79.
Fu KT, Ohba H. Crassulaceae. In: Wu CY, Raven PH, editors. Flora of China. Beijing: Science Press; 2001. p. 202–68.
Rohloff J. Volatiles from rhizomes of Rhodiola rosea L. Phytochemistry. 2002;59:655–61.
Svenning JC, Fløjgaard C, Marske KA, Nógues-Bravo D, Normand S. Applications of species distribution modeling to paleobiology. Quaternary Sci Rev. 2011;30:2930–47.
Alvarado-Serrano DF, Knowles LL. Ecological niche models in phylogeographic studies: applications, advances and precautions. Mol Ecol Resour. 2014;14:233–48.
Zhang JQ, Meng SY, Rao GY. Phylogeography of Rhodiola kirilowii (Crassulaceae): a story of miocene divergence and quaternary expansion. PLoS One. 2014;9:e112923.
Mayuzumi S, Ohba H. The phylogenetic position of eastern Asian Sedoideae (Crassulaceae) inferred from chloroplast and nuclear DNA sequences. Syst Bot. 2004;29:587–98.
Taberlet P, Gielly L, Pautou G, Bouvet J. Universal primers for amplification of 3 noncoding regions of chloroplast DNA. Plant Mol Biol. 1991;17:1105–9.
Hamilton MB. Four primer pairs for the amplification of chloroplast intergenic regions with intraspecific variation. Mol Ecol. 1999;8:521–3.
Liu PL, Wan Q, Guo YP, Yang J, Rao GY. Phylogeny of the genus Chrysanthemum L.: evidence from single-copy nuclear gene and chloroplast DNA sequences. PLoS One. 2012;7:e48970.
Thompson JD, Higgins DG, Gibson TJ. Clustal-W - improving the sensitivity of progressive multiple sequence alignment through sequence weighting, position-specific gap penalties and weight matrix choice. Nucleic Acids Res. 1994;22:4673–80.
Librado P, Rozas J. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics. 2009;25:1451–2.
Swofford DL. PAUP*. Phylogenetic analysis using parsimony (* and other methods). Version 4;2004.
Felsenstein J. Confidence-limits on phylogenies - an approach using the bootstrap. Evolution. 1985;39:783–91.
Posada D, Buckley TR. Model selection and model averaging in phylogenetics: advantages of akaike information criterion and Bayesian approaches over likelihood ratio tests. Syst Biol. 2004;53:793–808.
Posada D. ModelTest server: a web-based tool for the statistical selection of models of nucleotide substitution online. Nucleic Acids Res. 2006;34:W700–3.
Ronquist F, Huelsenbeck JP. MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics. 2003;19:1572–4.
Bandelt HJ, Forster P, Rohl A. Median-joining networks for inferring intraspecific phylogenies. Mol Biol Evol. 1999;16:37–48.
Drummond AJ, Rambaut A. BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol Biol. 2007;7:214.
Richardson JE, Pennington RT, Pennington TD, Hollingsworth PM. Rapid diversification of a species-rich genus of neotropical rain forest trees. Science. 2001;293:2242–5.
Yamane K, Yano K, Kawahara T. Pattern and rate of indel evolution inferred from whole chloroplast intergenic regions in sugarcane, maize and rice. DNA Res. 2006;13:197–204.
Zhang Q, Yang R, Wang Q, Liu JQ. Phylogeography of Juniperus przewalskii (Cupressaceae) inferred from the chloroplast DNA trnT-trnF sequence variation. Acta Phytotaxon Sin. 2005;43:503–12.
Excoffier L, Lischer HEL. Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and windows. Mol Ecol Resour. 2010;10:564–7.
Dupanloup I, Schneider S, Excoffier L. A simulated annealing approach to define the genetic structure of populations. Mol Ecol. 2002;11:2571–81.
Tajima F. Statistical-method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics. 1989;123:585–95.
Fu YX. Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics. 1997;147:915–25.
Rogers AR, Harpending H. Population-growth makes waves in the distribution of pairwise genetic differences. Mol Biol Evol. 1992;9:552–69.
Schneider S, Excoffier L. Estimation of past demographic parameters from the distribution of pairwise differences when the mutation rates very among sites: application to human mitochondrial DNA. Genetics. 1999;152:1079–89.
Rogers AR. Genetic-evidence for a pleistocene population explosion. Evolution. 1995;49:608–15.
Excoffier L, Laval G, Schneider S. Arlequin (version 3.0): an integrated software package for population genetics data analysis. Evol Bioinforma. 2005;1:47–50.
Harpending HC. Signature of ancient population-growth in a low-resolution mitochondrial-DNA mismatch distribution. Human Biol. 1994;66:591–600.
Wolfe KH, Li WH, Sharp PM. Rates of nucleotide substitution vary greatly among plant mitochondrial, chloroplast, and nuclear DNAs. Proc Natl Acad Sci U S A. 1987;84:9054–8.
Phillips SJ, Dudík M. Modeling of species distributions with Maxent: new extensions and a comprehensive evaluation. Ecography. 2008;31:161–75.
Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A. Very high resolution interpolated climate surfaces for global land areas. Int J Climatol. 2005;25:1965–78.
Peterson AT, Papeş M, Soberón J. Rethinking receiver operating characteristic analysis applications in ecological niche modeling. Ecol Model. 2008;213:63–72.
Elith J, Leathwick JR. Species distribution models: ecological explanation and prediction across space and time. Annu Rev Ecol Evol Syst. 2009;40:677–97.
Fielding AH, Bell JF. A review of methods for the assessment of prediction errors in conservation presence/absence models. Environ Conserv. 1997;24:38–49.
Hijmans RJ, Guarino L, Cruz M, Rojas E. Computer tools for spatial analysis of plant genetic resources data: 1. DIVA-GIS Plant Genetic Resources Newsletter. 2001;127:15–9.
Meng LH, Yang R, Abbott RJ, Miehe G, Hu TH, Liu JQ. Mitochondrial and chloroplast phylogeography of Picea crassifolia Kom. (Pinaceae) in the Qinghai-Tibetan plateau and adjacent highlands. Mol Ecol. 2007;16:4128–37.
Gao QB, Zhang DJ, Duan YZ, Zhang FQ, Li YH, Fu PC, Chen SL. Intraspecific divergences of Rhodiola alsia (Crassulaceae) based on plastid DNA and internal transcribed spacer fragments. Bot J Linn Soc. 2012;168:204–15.
Ren G, Mateo RG, Liu J, Suchan T, Alvarez N, Guisan A, Conti E, Salamin N. Genetic consequences of quaternary climatic oscillations in the Himalayas: Primula tibetica as a case study based on restriction site-associated DNA sequencing. New Phytol. 2017;213:1500–12.
Petit RJ, Aguinagalde I, de Beaulieu J-L, Bittkau C, Brewer S, Cheddadi R, Ennos R, Fineschi S, Grivet D, Lascoux M. Glacial refugia: hotspots but not melting pots of genetic diversity. Science. 2003;300:1563–5.
Avise JC. What is the field of biogeography, and where is it going? Taxon. 2004;53:893–8.
Zheng BX, Xu QQ, Shen YP. The relationship between climate change and quaternary glacial cycles on the Qinghai-Tibetan plateau: review and speculation. Quatern Int. 2002;97:93–101.
Dynesius M, Jansson R. Evolutionary consequences of changes in species’ geographical distributions driven by Milankovitch climate oscillations. PNAS. 2000;97:9115–20.
Jansson R, Dynesius M. The fate of clades in a world of recurrent climatic change: Milankovitch oscillations and evolution. Annu Rev Ecol Evol Syst. 2002;33:741–77.
Avise JC. Phylogeography: retrospect and prospect. J Biogeogr. 2009;36:3–15.
Liu JQ, Sun YS, Ge XJ, Gao LM, Qiu YX. Phylogeographic studies of plants in China: advances in the past and directions in the future. J Syst Evol. 2012;50:267–75.
Bennett KD. Continuing the debate on the role of quaternary environmental change for macroevolution. Philos Trans R Soc Lond Ser B Biol Sci. 2004;359:295–303.
Zhang ML, Uhink CH, Kadereit JW. Phylogeny and biogeography of Epimedium/Vancouveria (Berberidaceae): western north American - east Asian disjunctions, the origin of European mountain plant taxa, and east Asian species diversity. Syst Bot. 2007;32:81–92.
Brochmann C, Brysting AK. The Arctic – an evolutionary freezer? Plant Ecol Divers. 2008;1:181–95.
Levsen ND, Tiffin P, Olson MS. Pleistocene speciation in the genus Populus (Salicaceae). Syst Biol. 2012;61:401–12.
Wang Q, Abbott RJ, Yu QS, Lin K, Liu JQ. Pleistocene climate change and the origin of two desert plant species, Pugionium cornutum and Pugionium dolabratum (Brassicaceae), in Northwest China. New Phytol. 2013;199:277–87.
Pharmacopoeia Commission of the Ministry of Health of the People's Republic of China. Pharmacopoeia of the People’s Republic of China. Beijing: Chemical Industry Press; 2010.
Lei YD, Gao H, Tsering T, Shi SH, Zhong Y. Determination of genetic variation in Rhodiola crenulata from the Hengduan Mountains region, China using inter-simple sequence repeats. Genet Mol Biol. 2006;29:339–44.
Fordham DA, Brook BW, Moritz C, Nogues-Bravo D. Better forecasts of range dynamics using genetic data. Trends Ecol Evol. 2014;29:436–43.
Costion CM, Simpson L, Pert PL, Carlsen MM, Kress JW, Crayn D. Will tropical mountaintop plant species survive climate change? Identifying key knowledge gaps using species distribution modelling in Australia. Biol Conserv. 2015;191:322–30.
Pauli H, Gottfried M, Reiter K, Klettner C, Grabherr G. Signals of range expansions and contractions of vascular plants in the high Alps: observations (1994-2004) at the GLORIA master site Schrankogel, Tyrol, Austria. Glob Chang Biol. 2007;13:147–56.
The authors thank Yu-Qu Zhang for his assistance in sampling collecting, and Qian-Long Liang for his advice in ecological niche modeling. We appreciate the constructive suggestions from two anonymous reviewers, which made the paper much better.
This research was supported by National Natural Science Foundation of China (No. 31500177), the Fundamental Research Funds for the Central Universities (No. GK201603065), and Shaanxi Science and Technology Program (No. 2016JQ3025). The funding body didn’t participate any activities regarding the paper, including design of the study, sampling collection, analysis, data interpretation and in manuscript writing.
Availability of data and materials
The sequences used in this study have been deposited in the GenBank database (accession nos. MH371156 - MH371168). The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.
Jian-Qiang Zhang is an associate professor in Shaanxi Normal University who is interested in evolutionary history, adaptation and speciation of alpine plants in the Qinghai-Tibetan Plateau, especially in Crassulaceae.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Table S1. Bioclimatic variables (BIO1 to BIO19) from WorldClim . Variables marked with an asterisk (*) were used for the climatic niche models for Rhodiola crenulata. (DOCX 67 kb)
Table S2. Haplotype composition of 16 sampled populations of Rhodiola crenulata based on the cpDNA data set. (DOCX 16 kb)
Table S3. Haplotype composition of 16 sampled populations of Rhodiola crenulata based on the ITS data set. (DOCX 15 kb)
Figure S1. Divergence time of R. crenulata and its closest relative based on the plastid DNA haplotypes estimated with BEAST. Gray bars indicates 95% highest posterior density intervals. (TIF 93 kb)
Figure S2. Correlation between the F statistics and grouping number (K = 2–10) from the SAMOVA results. (a) results based on cpDNA haplotypes; (b) results based on ITS ribotypes. (TIF 887 kb)
Figure S3. Historical demography for overall populations and in each clade based on the plastid DNA dataset. Clade I and clade II correspond to the groups in the Bayesian phylogenetic tree in Fig. 4a. Mismatch distribution showing histogram of observed mismatch frequencies and best-fit curve of the sudden expansion model. (TIF 2056 kb)
About this article
Cite this article
Zhang, YZ., Zhu, RW., Zhong, DL. et al. Nunataks or massif de refuge? A phylogeographic study of Rhodiola crenulata (Crassulaceae) on the world’s highest sky islands. BMC Evol Biol 18, 154 (2018). https://doi.org/10.1186/s12862-018-1270-6
- Massif de refuge
- Qinghai-Tibetan plateau
- Quaternary climatic oscillations
- Quaternary speciation