Gone with the trees: Phylogeography of Rhodiola sect. Trifida (Crassulaceae) reveals multiple refugia on the Qinghai-Tibetan Plateau

Accepted Manuscript Gone with the trees: Phylogeography of Rhodiola sect. Trifida (Crassulaceae) reveals multiple refugia on the Qinghai-Tibetan Plateau Yuan-Cong Li, Da-Lv Zhong, Guang-Yuan Rao, Jun Wen, Yi Ren, Jian-Qiang Zhang PII: DOI: Reference:

S1055-7903(17)30602-4 https://doi.org/10.1016/j.ympev.2018.01.001 YMPEV 6020

To appear in:

Molecular Phylogenetics and Evolution

Received Date: Revised Date: Accepted Date:

17 August 2017 31 December 2017 2 January 2018

Please cite this article as: Li, Y-C., Zhong, D-L., Rao, G-Y., Wen, J., Ren, Y., Zhang, J-Q., Gone with the trees: Phylogeography of Rhodiola sect. Trifida (Crassulaceae) reveals multiple refugia on the Qinghai-Tibetan Plateau, Molecular Phylogenetics and Evolution (2018), doi: https://doi.org/10.1016/j.ympev.2018.01.001

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Gone with the trees: Phylogeography of Rhodiola sect. Trifida (Crassulaceae) reveals multiple refugia on the Qinghai-Tibetan Plateau

Yuan-Cong Lia, Da-Lv Zhonga, Guang-Yuan Raob, Jun Wenc, Yi Rena, Jian-Qiang Zhanga* a

College of Life Sciences, Shaanxi Normal University, Xi’an 710119, China

b

c

College of Life Sciences, Peking University, Beijing 100871, China

Department of Botany, National Museum of Natural History, MRC 166, Smithsonian

Institution, Washington, D.C. 20013-7012, U. S. A. *Corresponding author: [email protected] (J. Q. Zhang) Running title: Phylogeography of Rhodiola sect. Trifida

1

ABSTRACT Quaternary climatic oscillations have had tremendous effects on current distribution of species. Previous studies unraveled multiple microrefugia on the Qinghai-Tibetan Plateau (QTP) in two woody plants. Still we know little whether herbs growing in forests responded to climatic oscillations similarly. We herein conducted a phylogeographic study on Rhodiola sect. Trifida, an herbaceous group endemic to the QTP, which mainly growing on the forest floors, using plastid and ITS sequences as well as ecological niche modeling. The origin and divergence of major clades of sect. Trifida were in accordance with the last phase of the QTP uplifts. Mismatch distribution analysis indicated a range expansion dated to ca. 135 thousand years ago. A high frequency and an even distribution of private haplotypes in both plastid and ITS data sets throughout the distribution of sect. Trifida were detected. The ecological niche modeling results showed that there were suitable habitats on the QTP platform during the LGM. Our results found that multiple microrefugia existed on the QTP platform, supporting the hypothesis that species with similar geographic distribution and inhabiting the same community had similar responses to the Quaternary climatic oscillations. Furthermore, species delimitations in sect. Trifida need to be tested based on integrative evidence from morphological, ecological and genetic data. Keywords: Microrefugia; Phylogeography; Qinghai-Tibetan Plateau; Quaternary climatic oscillations; Rhodiola sect. Trifida; species delimitation.

2

1. Introduction Phylogeography deals with the spatial arrangement of genetic lineages within species or among closely related species (Avise, 2009). Accumulated evidence has suggested that the Quaternary climatic oscillations have had tremendous effects on current distribution of species (Hewitt, 2000; Hewitt, 2004). The effects of the Quaternary climatic oscillations on species ranges may have varied with latitudes and topography (Hewitt, 2004). Of great interest are the mountainous regions with varying altitudes, as montane species might have descended as the climate cooled and survived during the oscillations in the same region. Mountainous regions are also rich in biodiversity with the orogenic movements, complex topography and existence of complex microhabitats (Hewitt, 2004; Hughes and Atchison, 2015). Phylogeographic studies on these mountain areas may provide important insight into how geology and climatic oscillations have shaped the current biodiversity and distributional pattern of alpine plants. As the earth’s largest plateau, the Qinghai-Tibetan Plateau (QTP) was formed by several uplift events after the collision of the Indian plate with the Eurasian plate about 40 Mya (million years ago) (Harris, 2006; Guo et al., 2002; Spicer et al., 2003). Harboring more than 12,000 seed plant species, the QTP and its adjacent areas (the “QTP” hereafter means both the QTP platform and its adjacent areas) have been a hotspot for many phylogenetic and phylogeographic studies (Wen et al., 2014). Previous studies have revealed general patterns of plants on the QTP responding to the 3

glacier climatic oscillations: contracting to Hengduan Mountains during glaciations and recolonizing onto the plateau (Yang et al., 2008; Zhang et al., 2010); surviving in refugia that existed on the plateau platform during glaciations and local expansion during interglacial or postglacial periods (Wang et al., 2009a, b); surviving in multiple microrefugia during the Last Glacial Maximum (LGM, ca. 0.02 Mya) throughout the current distribution of species, as confirmed by two previous studies in woody plants on the QTP (Opgenoort et al., 2010; Wang et al., 2010); or as suggested by Luo et al. (2016) that geographic distribution show little or no impact of the Quaternary glaciations in some species. Although these studies demonstrated some important patterns, they have mainly focused on woody plants or herbs on alpine meadows. Still we know little how herbaceous species growing in forests responded to the glacial climatic oscillations, as they should have been affected greatly by the range shifts of forests. Most previous phylogeographic studies only sampled a single species at the population level. This sampling strategy may show severe bias in phylogeographic inference because introgressive gene flow, which results in chloroplast capture, may be common in closely related species (Rieseberg and Soltis, 1991; Yi et al., 2015). A haplotype revealed in a species or a population, whether it is rare or common, might be newly evolved or inherited from a common ancestor; and it also could be obtained by interspecific gene flow from other closely related species. This kind of bias should be taken into account more carefully in most plastid DNA marker based phylogeographic 4

studies. Several studies drew different conclusions after more closely related species were sampled compared to previous single-species based studies (e.g., Qin et al., 2013; Su et al., 2010). Our study focuses on Rhodiola sect. Trifida, which is endemic to the QTP and its adjacent area, and is composed of 3-9 species, according to different authors. Species of this section are characterized by bisexual white flowers, short flowering stems, and coarsely lobulated leaves with relatively large rounded teeth in Rhodiola (Fu and Fu, 1984; Ohba, 1978). A molecular phylogenetic study based on plastid and nuclear ribosomal internal transcribed spacer (ITS) sequences supported its monophyly (Zhang et al., 2014a, b). Ohba (1978) only included three species of this section, i.e., R. chrysanthemfolia, R. sinuata and R. liciae, while Fu and Fu (1984) delimited nine species in this section mainly based on leaf size and degree of leaf division (see Table S1). Yet these leaf characters could be easily affected by environmental factors. The recent Flora of China treatment largely followed Fu and Fu (1984), and included eight species in the section (Fu and Ohba, 2001). Species of sect. Trifida have a center of distribution in the southern QTP and the Hengduan Mountain region, with some populations of R. sacra and R. alterna scattered across the interior of the QTP. Species of this section mainly grow on the forest floors, in alpine shrublands, at forest margins, and even on tree trunks. Rhodiola sacra and R. sinuata can expand to rock crevices on grassland slopes. Previous phylogeographic studies on juniper forests revealed that a highest tree line existed at LGM in the southern Tibetan Plateau (Opgenoorth et al., 5

2010). Thus, it is reasonable to hypothesize that there are also multiple microrefugia that existed in LGM for species of sect. Trifida, as the juniper forests may have provided suitable habitats for them to survive during the LGM. Here we conducted a phylogeographic study of Rhodiola sect. Trifida employing the maternally inherited cpDNA intergenic spaces trnL-F, trnS-G, and psbA-trnH, as well as biparentally inherited ITS sequences. We also used ecological niche modeling (ENM) to examine the potential range shifts of sect. Trifida in response to glacial climatic oscillations. Through the above-mentioned analyses, we aim to test the hypothesis that multiple micro-refugia on the southern part of the QTP existed for species of this section to survive during LGM, and to infer the historical geological and climatic factors that may have led to the current distribution and genetic variation pattern of the QTP and its adjacent areas endemic lineage Rhodiola sect. Trifida. 2. Materials and methods 2.1 Population sampling Fresh leaves were collected throughout the distribution range of all species in Rhodiola sect. Trifida from 2009 to 2016. We collected a total of 590 individuals from 43 populations, representing most taxa in previous taxonomic treatments (except Rhodiola tieghemii (Hamet) S. H. Fu and R. dielsiana (Limpr. f.) S. H. Fu, which should be merged to R. chrysanthemifolia (H. Lév.) S. H. Fu according to Fu and Ohba (2001)) (Figs. 1, 6D; Table S1). Sample size per population ranged from 2 to 20 depending on population size. Leaves from each individual were dried immediately in 6

silica gel and stored at room temperature. Each individual collected was at least 15 m apart to avoid sampling bias toward closely related individuals. 2.2 DNA extraction, PCR amplification and sequencing Total genomic DNA was extracted from 30 mg of silica-dried leaves using a modified CTAB method (Doyle, 1991). Three cpDNA intergenic spacers and the nrDNA ITS fragment were amplified to assess phylogeographic patterns using polymerase chain reaction (PCR). Primers used were psbAF (5’-GTTATGCATGAACGTAATGCTC- 3’) and trnHR (5’-CGCGCATGGTGGAT TCAC AAATC-3’) for psbA-trnH (Hamilton, 1999), trnG (GAACGAATCACACTTTTACCAC) and trnS (GCCGCTTT AGTCCA CTCAGC) for trnS-G (Hamilton, 1999), c (5’-CGAAATC GGTAGA CGCTACG -3’) and f (5’-ATTTGAACTGGTGACACGAG-3’) for trnL-F (Pierre et al., 1991), and ITS1 (5’-TCC GTAGGTGAACCTGCGG- 3’) and ITS4 (5’-TCCTCC GCTTATTG ATATGC-3’) for ITS (White et al., 1990). Each PCR amplification was implemented in a 20 μl volume containing 10 μl of 2 × Taq master mix (blue) reaction buffer, 0.8 μl of each primer at 100 ng/μl, 1 μl of genomic DNA at 50 ng/μl, and 7.8 μl of Milli-Q H2O. The cycling profile of the ITS region was: initial denaturation at 95°C for 4 minutes (min), followed by 38 cycle of 30 seconds at 95°C, 30 seconds of annealing at 55°C, 30 seconds of elongation at 72°C, and ending with a 7 min extension at 72°C. The annealing temperature of amplification set for each pair of the plastid primers ranged from 50°C to 52°C. PCR fragments were visualized on agarose gels and purified using 7

polyethylene glycol (PEG) precipitation prior to sequencing. Most of the ITS samples were directly sequenced, although a few individuals (21 out of 493) had ambivalent peaks in the chromatogram. These individuals were ligated onto pGEM-T Easy Vector using a Promega Kit (Promega Corporation, Madison, WI, USA). Plasmids containing the right insertion were chosen for sequencing reactions. Sequencing was carried out by Sangon Biotech (Shanghai) Co., Ltd., China, using an ABI 3730XL sequencer (Applied Biosystems, Foster City, CA, USA). Primers used for sequencing were the same as those for amplifications. Contigs were assembled and edited using Geneious R9 (Masters et al., 2011) (http://www.geneious. com/). The sequences were deposited in GenBank under accession numbers MG309161-MG309305. Sequences were aligned using MUSCLE 3.8.31 (Edgar, 2004). The three plastid fragments (psbA-trnH, trnL-F and trnS-G) were concatenated for the analysis. 2.3 Population genetics For each sequence dataset, haplotype diversities (h) and nucleotide diversities (π) were calculated using ARLEQUIN 3.5 (Excoffier et al., 2005). These measurements of molecular diversity reflect both the overall number of haplotypes and their molecular divergence (Guest et al., 2014). Differentiation within populations, among populations within groups and among groups, was calculated by analyses of molecular variance (AMOVA) to estimate genetic variation in different species implemented in ARLEQUIN 3.5, with significance tested via 10,000 permutations (Excoffier and 8

Heckel, 2006). The highest ranking (among groups) was divided by species delimitation. To test if there is phylogeographic structure, we calculated four measures of genetic differentiation: divergence between populations GST and NST, average gene diversity with a population (HS), and total gene diversity (HT). GST is a differentiation measure that is based on allele frequencies only, while NST takes into account the similarities between haplotypes (i.e., the number of mutations between haplotypes). Thus, a higher NST than GST indicates a phylogeographic structure with closely related haplotypes being more likely to co-occur close to each other. The two parameters were compared using a permutation test with 1000 permutations and the U-statistic as implemented in PERMUT (http://www. Pierroton.inra.fr/genetics/labo/Software/Permut). To identify signatures of demographic history, the mismatch distribution and neutrality tests (Tajima’s D and Fu’s Fs) of each clade recovered based on the plastid dataset were performed using ARLEQUIN 3.5. If an expansion hypothesis was true, significant negative values of Tajima’s D and Fu’s FS statistics could be observed. We used all haplotypes for each clade to conduct the mismatch distribution analysis. We used 1000 parametric bootstrap replicates to test the fitness between observed mismatch distributions and the theoretical distribution under a recent spatial expansion model, examining the sum of squared deviations (SSD) for observed and expected mismatch distributions, and the raggedness index (HRag) of Harpending (Harpending, 1994). If a group was identified to have experienced expansion, we used the parameter-value for the mode of the 9

mismatch distribution (τ) to estimate time (in generations) since expansion with the equation t = τ/2u (Rogers and Harpending, 1992) (Rogers, 1995). Here, u was calculated as (u = μ × k × g), where μ is the number of substitutions per site per year (s s-1 yr-1), k is the average sequence length in the present study and g is the generation time in years (age when first reproduction occurs). In the present study, our value for k was 1801 bp. Chloroplast DNA substitution rates for most angiosperm species have been estimated to be in the range 1-3×10-9 substitutions per site per year (s s-1 yr-1) (Wolfe et al., 1987). As the two fragments we used are non-coding regions of the cpDNA genome, we assumed an evolutionary rate of 1.52×10 -9 s s-1 yr-1 for the plastid dataset (Yamane et al., 2006). Ten years were used as an approximation for g (J. Q. Zhang, pers. obs.). To test whether there is significance of isolation by distance between populations, the Mantel test with 1, 000 random permutations on matrices of pairwise population FST values and the geographic distances was performed using ARLEQUIN ver. 3.5 (Excoffier and Lischer, 2010). 2.4 Phylogenetic relationships and lineage divergences We selected two species from Rhodiola sect. Prainia S. H. Fu (Fu and Ohba, 2001) as the outgroups, based on previous phylogenetic results (Zhang et al., 2014a). We utilized DnaSP v. 5.0 (Librado and Rozas, 2009) to estimate chloroplast haplotypes and ITS ribotypes, which were determined from both indels and nucleotide substitutions, respectively. Phylogenetic relationships were reconstructed based on maximum parsimony (MP) and Bayesian inference (BI). Insertion/deletion events (indels) were 10

coded according to Simmons and Ochoterena’s (Simmons and Ochoterena, 2000) ‘simple’ indel coding using the software SeqState (Muller, 2005). PAUP 4.0 (Swofford, 2002) was used to conduct the MP analyses, and all characters and transitions/transversions were considered as equally weighted. MULTREES and TBR branch swapping were performed for heuristic searches. Starting trees were created by random addition sequences and repeated 1000 times. Support for individual nodes was evaluated by bootstrapping (Felsenstein, 1985). The bootstrap analyses used 1000 replicates with ten replicates of random addition sequences and NNI branch swapping, and then a maximum of 1000 trees of the random addition replicates were saved. jModelTest v. 2.1.4 (Darriba et al., 2012) was used to choose substitution model by the AIC (Darriba et al., 2012; Guindon and Gascuel, 2003). The Bayesian inference was performed using MrBayes version 3.2.1 (nst = 6, rates = gamma for the plastid dataset; nst = 1 for the ITS dataset) (Ronquist and Huelsenbeck, 2003). The analyses started with a random tree and used four Metropolis-coupled Markov chain Monte Carlo chains (each with 5,000,000 generations), saving one tree every 1,000 generations. The first 25% trees were discarded as burn-in. We constructed cpDNA and ITS haplotype networks with NETWORK v4.2.0.1 (Bandelt et al., 1999) in order to detect genealogical relationships among sequences with shallow genetic divergences. We assumed both site mutations and indels have evolved with equal likelihood and each indel has originated independently of all other indels. 11

BEAST v2.0.1 (Drummond et al., 2012) was used to infer the divergence times between lineages. We used both ITS and plastid datasets to conduct dating analysis with an uncorrelated lognormal relaxed molecular clock model in it. We used the program MEGA v6.06 (Tamura et al., 1994) to insure the best-fit molecular clock hypothesis by comparing log likelihood values with and without enforced molecular clock. Divergence times and confidence intervals were then estimated (Drummond and Rambaut, 2007; Rambaut and Drummond, 2007), with four gamma categories and the birth–death process of tree prior. Because there were no fossil records reported of sect. Trifida (Rhodiola) or their close relatives, we used substitution rates to estimate divergence times. The substitution rate of nrITS in shrubs and herbs varies from 3.46×10-9 to 8.69×10-9 s s-1 yr-1 (Richardson et al., 2001). We thus used a normal distribution prior to cover these ranges with a 95% confidence interval because of uncertainty of the rates, which is a mean of 6.075×10-9, and SD of 1.590×10-9. The substitution rate for the plastid dataset was the same as in the mismatch distribution analysis. FigTree v. 1.3.1 (Rambaut, 2009) was used to display a tree with ages for each node and their 95 % highest posterior density (HPD). 2.5 Ecological niche modeling We used ENM to examine the potential range shifts of sect. Trifida in response to glacial climatic oscillations. Localities of sect. Trifida were obtained from field collections by the authors and from online herbarium records (e.g., Chinese Virtual Herbarium and Global Biodiversity Information Facility). Data initially obtained from 12

online databases were double-checked with actual specimens to exclude misidentifications. MAXENT v. 3.3.3 (Phillips and Dudík, 2008) was used to model the ecological niches of species in this section. To ensure that there are no autocorrelations among the locations we used, we firstly removed sites that are too close to each other manually. Also, we checked “remove duplicate presence records” in the MAXENT program. A total of 203 spatially unique localities were used for the analysis (Table S5). Environmental layers of 19 bioclimatic variables (Table S4) for the Last Glacial Maximum (LGM, ca. 20 kya, kilo years ago) (MIROC: Hasumi and Emori, 2004; CCSM: Collins et al. 2006) and the present at a spatial resolution of 2.5 arc-minutes, as well as Last Interglacial (LIG, ca. 130 kya) at a resolution of 30 arc-second (Otto-Bliesner et al. 2006) were downloaded from the WorldClim website and employed for the modeling (Hijmans et al., 2005). To exclude highly correlated climate variables, pairwise correlations were examined among the 19 variables within the distributional area of sect. Trifida. The six variables (Table S4) with pairwise Pearson correlation coefficients below 0.7 were used. Area under the “receiver operating characteristic (ROC) curve” (AUC) (Elith and Leathwick, 2009; Peterson et al., 2008) values were used to evaluate the accuracy of each model prediction. The threshold for good performance was set to 0.7 (Fielding and Bell, 1997). DIVA-GIS version 7.5 (Hijimans et al., 2001) was used to draw the range of suitable distributions, and to resample the LIG bioclimatic layers to achieve the same resolution with LGM and present. 13

3. Results 3.1 Plastid DNA sequence variation and haplotype distribution The length of aligned plastid sequences was 325 bp, 621 bp, 855 bp for psbA-trnH, trnS-G, and trnL-F, respectively. The concatenated plastid dataset had the alignment length of 1801 bp. Across all 590 sampled individuals, 55 haplotypes were identified (H1-H55), characterized by 61 variable sites, of which 58 were parsimony informative, and two indels (9 bp and 23 bp, respectively). The number of haplotypes per population ranged from one to six (Fig. 1A; Table S1). The two most frequent haplotypes (H2, H25) accounted for 29% of the samples. H2 was the most common haplotype, which occurred in 16 of the 43 populations (Fig. 1A). The majority of haplotypes (40 out of 55) were fixed in one population (private haplotypes). Strikingly, these private haplotypes were evenly spread all over the distributional range (Fig. 1B). At the sectional level, the estimated haplotype diversity was h = 0.923. Haplotype diversity varied across populations, ranging from 0.000 to 0.929, with the population in Lazi county of R. sacra (Sa_LZX) having the highest h value (Table S1). Nucleotide diversity π was 0.00346 at the sectional level and ranged from 0.000 to 0.0046 in different populations, with population Sa_LZX having the highest π value. Nearly all populations from southern Tibet exhibit above-average cpDNA genetic diversity (Table S1). Within-population gene diversity HS (0.314) was much lower than the total gene diversity HT (0.936). Of the 43 sampled populations, 24 harbored (55.8%) private haplotype(s) (Fig. 1B). 14

3.2 Nuclear ITS sequence variation and ribotype distribution In the 21 individuals that had ambivalent peaks in the chromatogram, we confirmed via cloning sequencing that 9 of them (Sa_GBJ_11; Sa_MK_3_9; Sa_NML_14; Ch_SKXS_1; Ch_LZ_19; Al_YA_2; Ov_MXG_1_1; Ov_MXG_1_4; Se_QSX_1_8; See details in Table S1) had two ribotypes. A total of 47 (R1–R47) different ITS sequences (ribotypes) were detected from the 493 individuals, determined by 110 variable sites. The length of aligned ITS sequences was 648 bp, with a length variation from 633–648 due to low sequence quality in both ends in some individuals. 27 of the 47 ribotypes occurred only in a single population, while the other 20 were shared by at least two populations (Fig. 2). Similar to cpDNA haplotypes, these private ribotypes were also evenly spread throughout the distributional range (Fig. 2B). The most widespread ribotype was R7, which was found in eight populations (Fig. 2A). Ribotype diversity was estimated to be h = 0.948 at the sectional level, ranging from 0.000 to 0.900 in different populations. Population Sa_LZX_1 in Lazi County has the highest h value. At the sectional level, nucleotide diversity was 0.023, ranging from 0 to 0.032 in different populations. Populations from southern Tibet also exhibit above average ITS genetic diversity (Table S1). In the ITS datasets, the value of total genetic diversity (HT = 0.959) was significantly higher than the average within-population diversity (HS = 0.419). 3.3 Phylogenetic relationships and lineage divergences

15

For both the plastid and the ITS datasets, our phylogenetic analyses based on MP and Bayesian methods yielded largely congruent topologies (Fig. 3A; Fig. 4A). The monophyly of sect. Trifida was well supported in both datasets. For the cpDNA dataset, we found two well-supported groups, with the remaining haplotypes forming a polytomy (Fig. 3A). Clade I included haplotypes H11, H21-H23, H31, H37, H47 and H51. Clade II comprised haplotypes H5–H7, H17–H18, H20, H32–H33 and H55. Clade I was distributed in the Hengduan Mountains, and Clade II was in both the Hengduan Mountains and the QTP platform (Fig. 1A). The remaining haplotypes had a wide distributional range across the Hengduan Mountains and the QTP platform. The network diagram was consistent with the phylogenetic relationships (Fig. 3B), but showed relationships between haplotypes in more detail. For the ITS ribotypes, we detected four clades (Fig. 4A): the main clade A was divided into subclade A1 and subclade A2 (Fig. 4A). The two ribotypes of subclade C were private, while other subclades had both private and shared ribotypes. The ribotype network (Fig. 4B) depicted the relationships among the ribotypes in more detail. We dated the divergence times of sect. Trifida using the ITS dataset, which suggested that the separation between sect. Trifida and the outgroups was at 3.54 Mya (Node 0; 95%HPD: 2.39–4.58 Mya), and the divergence time for the two major clades was at 3.12 Mya (Node 1; 95%HPD: 2.22–4.10 Mya) in the late Pliocene (Fig. 5). Further diversification of the multiple main ribotypes (subclades) took place at around

16

1.5–1.7 Mya during the Pleistocene (Fig. 5). The dating analysis based on the plastid DNA markers yielded congruent results with the ITS dataset (Fig. S1). 3.4 Population structure Both cpDNA and ITS datasets revealed that the values of total genetic diversity (HT) were significantly higher than those of the average within-population diversity (HS) (cpDNA, HT = 0.936, HS = 0.314; ITS, HT = 0.959, HS = 0.419). For the plastid dataset, a strong phylogeographic structure was detected, as the NST (0.853) was significantly larger than the GST (0.665, p < 0.01), suggesting that gene flow is lower than mutation rate. Nevertheless, NST (0.637) was not significantly larger than the GST (0.564) for the ITS dataset. The difference between ITS and cpDNA was also found in each regional group (Table 1). Analyses of molecular variance (AMOVA) showed that 19.30% of the total plastid DNA genetic variation was found among groups, and 66.51% among population within groups, and 14.19% within population (Table 2). For the ITS data, 8.63% of the total genetic variation was found among groups, and 60.91% among populations within groups, and 30.48% within the population (Table 2). The Mantel test based on cpDNA and ITS suggested a non-significance pattern of isolation-by distance (cpDNA: r = 0.1407, p > 0.05; ITS: r = -0.121, p > 0.05). 3.5 Demographic analyses The mismatch distribution analysis for the cpDNA dataset revealed that the mismatch distribution for all samples was unimodal, while clade I and clade II showed a multimodal mismatch distribution (Fig. S2). We also found that the observed variance 17

(SSD = 0.005, p > 0.05) and the raggedness index (HRag = 0.106, p > 0.05) for all samples were not significantly different from the expectation under the population expansion model (Table 3). Tajima’s D (-0.848, p > 0.05) and Fu’s Fs ( -12.709, p < 0.05) for all samples were significantly negative, while those of clade I and clade II present no significant negative values. Thus, all evidence strongly indicated that species of sect. Trifida underwent a recent spatial expansion and the estimated time of this expansion was around 0.135 (0.068-0.205) Mya (Table 3). 3.6 Ecological niche modeling Ecological niche modeling performed well in prediction of geographic distribution of species in sect. Trifida, with the area under the receiver operating characteristic curve = 0.983 ± 0.005 (mean ± SD). This means that our models is distinct from random expectation. As expected, the distributions of species in sect. Trifida during the LIG (ca. 130 kya) were smaller than that during the LGM (ca. 20 kya), almost confined to the Hengduan Mountain area (Fig. 6A, B, C). Both LGM models (CCSM and MIRCO) generated potential distributions that are largely consistent with each other (Figs. 6B, C). During the LGM, both the CCSM (Fig. 6B) and the MIROC (Fig. 6C) models showed a distribution pattern more or less the same as that in the present day, but with smaller distribution area on the plateau platform. Notably, the data showed that there were suitable habitats during the LGM in the Himalayan area (Fig. 6B, C). Furthermore, the current distributional predictions were largely congruent with the actual distribution of the species (Fig. 6D). 18

4. Discussion 4.1 Microrefugia for Rhodiola sect. Trifida during glaciations In the last two decades, species responding to the Quaternary climatic oscillations on the QTP and its adjacent areas have been extensively studied (Gao et al., 2012; Meng et al., 2007; Opgenoorth et al., 2010; Ren et al., 2017; Wang, 2009a, b; Wang et al., 2010; Yang et al., 2008). A recent study on four subnival species of Himalaya-Hengduan Mountains showed that Quaternary glaciations may have little effect on the distribution of studied species (Luo et al., 2016). On the other hand, Pedicularis longiflora (Yang et al., 2008) and Picea crassifolia (Meng et al., 2007) demonstrated declined genetic diversity or even single haplotypes on the QTP platform but a high level of genetic diversity and private haplotypes in the Hengduan Mountains. These data are thus consistent with the tabula rasa hypothesis (Nordal, 1987). Some other species like Potentilla glabra (Wang,et al., 2009a) and Aconitum gymnandrum (Wang, et al., 2009b) showed a different pattern that refugia existed both on the platform and the Hengduan Mountains. More strikingly, another pattern was exemplified by the Juniperus tibetica complex (Opgenoorth et al., 2010) and Hippophae tibetana (Wang et al., 2010), with an even distribution of populations with high genetic diversity and frequency of private haplotypes across the current distributional range, indicating “microrefugia” existed across the distributional range. The distribution of the Juniperus tibetica complex, which included five tree species of Juniperus, as well as the shrub species Hippophae tibetana is largely congruent with 19

that of sect. Trifida. Multiple LGM microrefugia throughout the current distributional range of the J. tibetica complex, partly well above 3500 m, marked the world’s highest known tree line. The same pattern as revealed in Hippophae tibetana demonstrated even higher LGM microrefugia (above 4000 m). As species of sect. Trifida often inhabit the forest floors, alpine shrublands, forest margins, and even tree trunks, they might have also endured glaciations on the QTP as the above-mentioned woody taxa might have provided them with the potential microhabitats. Our results were in accordance with the microrefugia hypothesis exemplified by the Juniperus tibetica complex and Hippophae tibetana. Firstly, 40 of the 55 plastid haplotypes and 27 of the 47 ribotypes were fixed in one population (Table S1; Figs. 1 and 2). These private haplotypes were evenly spread all over the distributional range (Figs. 1B, 2B), which were highly unlikely the results of dispersals from one macro-refugium. They could have evolved in situ in a minimum time span of 0.37 Mya, in the case of only one mutation between different private haplotypes and the substitution rate of 1.52×10-9 substitutions per site per year for cpDNA (Yamane et al., 2006). Our dating analysis based on the cpDNA dataset also showed that the divergence time between closest private haplotypes was at around 0.25 Mya (Fig. S1), which is well before the LGM. Besides, overall degrees of diversity are not significantly different between the plateau populations and the Himalayan populations or the Hengduan Mountain populations. Nearly all populations from southern Tibet exhibit above-average cpDNA genetic diversity (Table S1). This scenario is further 20

supported by our ENM results (Fig. 6): during LGM, there were still suitable habitats across the Himalayan region in both the CCSM and MIROC models (Fig. 6B, C). However, the ENM results should be treated with caution, as in mountainous areas, the 2.5 arc-minutes resolved bioclim data could be problematic because of the complex topography, especially in QTP, where plants could be easily affected by micro-habitats. All these results strongly support the microrefugia hypothesis for sect. Trifida. Thus, our results confirmed the long held hypothesis by comparative phylogeographers (e.g., Avise, 2009) that species with similar geographic distributions and inhabiting in the same community might have responded similarly to the Quaternary climatic oscillations. 4.2 Heterogeneous contributions of geography and climatic oscillations to the genetic structure of sect. Trifida As the world’s largest and highest plateau, the southern and southeastern parts of the QTP (the Himalaya and the Hengduan Mountains) were listed as the Earth’s biodiversity hotspots with numerous endemic species (Mittermeier et al., 2001, 2005; Myers et al., 2000). Species endemic to the QTP its adjacent areas may have been affected by both historical sequential uplifts of the mountains and the Quaternary climatic oscillations (Wen et al., 2014). In the present study, all species of sect. Trifida are endemic to the QTP and its adjacent areas (Ohba, 1987). In our BEAST analysis based on the ITS dataset, the origin of this section was dated to 3.54 Mya (95%HPD: 2.39–4.58 Mya), which was in accordance with the last phase of the QTP uplift (Li et 21

al., 2001). This result is consistent with that of a previous generic-level study, which hypothesized that the origin and diversification of Rhodiola were triggered by the intense uplift of the QTP (Zhang et al., 2014a). However, the relationship between uplift of the QTP and divergence of plant groups, especially young lineages, should be treated with caution. As Renner (2016) stated, there are evidence that the QTP got its height well before the Pliocene. Nevertheless, our estimation coincided with the range of most sectional level studies of the QTP region (summarized by Renner, 2016). Previous studies also demonstrated that both geological factors and climatic oscillations have deeply affected the distribution of biota of the QTP in both animals (Zhou et al., 2017) and plants (Li et al., 2013; Liu et al., 2014). In the case of sect. Trifida, geological factors, i.e., the intensive uplift of the QTP, could have acted as drivers of divergence of the major clades of the ITS tree, which was dated as around 1.4-3.1 Mya (Fig. 5). The heterogeneity of topography, climate and ecology on the QTP would have provided opportunities for species divergence and speciation (also see Zuo et al., 2015). Ribotypes of clades B and C are largely confined to the plateau platform, except R5 of population 43 (Ch_SKXS). However, ribotypes of clades A and D were scattered across the distributional area of sect. Trifida (Fig. 2). This pattern might be caused by extensive gene flow by pollen transfer after their divergence. Our Mantel test results based on the ITS sequences suggested no significance of the isolation-by-distance pattern (r = -0.121, p > 0.05), which confirmed strong pollen

22

mediated gene flow between populations. Besides, NST was not significantly larger than the GST for the ITS dataset, also suggesting the lack of a phylogeographic structure. For the plastid dataset, haplotypes of Clade I were restricted to the Hengduan Mountain area, and those of Clade II were in the Hengduan Mountains and on the QTP platform. Divergence time estimation suggested the divergence of the two groups as well as the other polytomy clades at ca. 2–4 Mya (Figs. S1), which was also concurrent with the last phase of the QTP uplift. It is noted that only H2 and H25 were widespread throughout the distributional area, and most of the haplotypes were restricted in a single population. This pattern indicated that gene flow through seeds was limited or even absent across the entire distributional range of species in this section. A strong phylogeographic structure was detected based on the plastid dataset, as the NST (0.853) was significantly larger than the GST (0.665, P < 0.01), thus conforming the limited gene flow via seeds. Furthermore, the AMOVA results showed that 66.51 % of the total plastid DNA variation was found among populations. In the mismatch distribution analysis, we did not detect signatures of range expansion in Clade I and Clade II, whereas a significant expansion was detected for all samples as a whole. Using 1.52×10-9 s s-1 yr-1 for the plastid dataset (Yamane et al., 2006), the estimated time of this expansion was at ca. 0.135 (0.068-0.205) Mya, well before the LGM and just at the beginning of the Penultimate Glaciation (Zheng et al., 2002). This expansion may have resulted in the widespread distribution of cpDNA haplotype H2 and H25. Our data thus provided further evidence that distributional pattern of species endemic to the QTP area 23

has been affected by both complex geological heterogeneity and climatic oscillations in the Quaternary. 4.3 Implications for species delimitations and phylogeography of close related species Species of sect. Trifida have been identified as constituting a monophyletic group in previous phylogenetic studies (Zhang et al., 2014a, b). The monophyly of the section is further supported by the following morphological synapormorphies: short flowering stems, and coarsely lobulated leaves with relatively large rounded teeth (Fu and Fu, 1984; Ohba, 1978). Nevertheless, there have been disagreements on the number of species in the section. Ohba (1978) recognized only three species, i.e., R. chrysanthemfolia, R. sinuata and R. liciae, while Fu and Fu (1984) delimited nine species in this section mainly based on leaf size and morphology, and degree of leaf division. Fu and Ohba (2001) subsequently recognized eight species, largely following Fu and Fu (1984). In Fu and Fu (1984)’s treatment, leaf distribution on the flowering stem, and the size and the margin of leaves are considered important diagnostic characteristics in delimitate other species of this section. However, these characters seem to be variable in our analysis. For example, the most pervasive haplotype H2 in the plastid data was shared by five species (Fig. 3B). The same pattern was also revealed in the ITS data (Fig. 4B). Another example comes from R. alterna, which was documented as distinct from its close relative R. sacra in leaf margin. The former had 2 shallowly incised serrations, while the latter had 4 or 5 dentate lobes on each side. In both datasets, R. alterna shared haplotypes with R. sacra (H1 in the plastid data and R7, 24

R24 and R25 in the ITS data; Figs. 3B and 4B). Thus, the two species may best be merged into one. Rhodiola liciae is distinct from other species of this section in its awn-like petal apex with a mucro (Fu and Ohba, 2001). In the plastid dataset, individuals of this species all possessed the private haplotype (H11), which showed 9 mutations apart from the nearest relatives (Fig. 3B), indicating a long-term separation from other species. However, the ITS network indicated pollen mediated gene flow between R. sacra and R. liciae (Fig. 4B). Furthermore, R. liceae is geographically confined to the Xishan Mountain of Kunming, Yunnan province. Based on all the above-mentioned evidence, the species status of R. liceae should be retained. It has been shown that the percentage of leaves with entire margins in living floras is linearly correlated with the average annual temperature, which in turn is affected by altitudes (Spicer et al., 2003). The degree of leaf division may be easily affected by environmental factors, and thus should not be used as diagnostic characteristics for species delimitations. Our data showed strong discrepancy with current species concept, thus further studies integrating morphology, ecology and genetic data needed to achieve a better species delimitation in this section (also see Spooner, 2016; Wen et al., 2015). Our study also demonstrated the importance of sampling more inclusive species in studying phylogeography of complex taxonomic groups. Species of these groups, like sect. Trifida in this case, are often closely related and mutually interfertile. Studying a single species will miss a lot of information regarding the whole genetic history of species group. 25

What’s more, closely related and interfertile species often have common genetic patterns reflecting closely shared phylogeographic histories. Thus we suggest that future phylogeographic studies sample more close related species to get a comprehensive picture of glacial history.

26

Acknowledgements We are grateful to Dr. Yu-Qu Zhang for his assistance of sample collection, and Dr. Gang Liu for ENM analysis. This study was supported by National Natural Science Foundation of China (no. 31500177), Shaanxi Science and Technology Program (no. 2016JQ3025), the Fundamental Research Funds for the Central Universities (no. GK201603065 to J.Q. Zhang and 2016CSZ006 to Y.C. Li).

27

References Avise, J.C., 2009. Phylogeography: retrospect and prospect. J. Biogeogr. 36, 3–15. Bandelt, H.J., Forster, P., Röhl, A., 1999. Median-joining networks for inferring intraspecific phylogenies. Mol. Biol. Evol. 16, 37–48. Collins, W.D., Bitz, C.M., Blackmon, M.L., Bonan, G.B., Bretherton, C. S., Carton, J.A., Chang, P., Doney, S., Hack, J.J., Kiehl, J.T., Henderson, T., Large, W.G., McKenna, D., Santer, B.D., Smith, R.D., 2006. The Community Climate System Model version 3 (CCSM3). J. Clim. 19, 2122–2143. Darriba, D., Taboada, G.L., Doallo, R., Posada, D., 2012. jModelTest 2: more models, new heuristics and parallel computing. Nat. Methods 9, 772–772. Doyle, J., 1991. DNA protocols for plants-CTAB total DNA isolation. In: Hewitt, G.M., Johnston, A.W., Young G.P.W. (Eds.), Molecular Techniques in Taxonomy. Springer Science & Business Media, Dordrecht, pp. 283–293. Drummond, A.J., Rambaut, A., 2007. BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol. Biol. 7, 214. Drummond, A.J., Suchard, M.A., Xie, D., Rambaut, A., 2012. Bayesian phylogenetics with BEAUti and the BEAST 1.7. Mol. Biol. Evol. 29, 1969–1973. Edgar, R.C., 2004. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic. Acids. Res., 32, 1792–1797.

28

Elith, J., Leathwick, J.R., 2009. Species distribution models: ecological explanation and prediction across space and time. Annu. Rev. Ecol. Evol. Syst. 40, 677–697. Excoffier, L., Heckel, G., 2006. Computer programs for population genetics data analysis: a survival guide. Nat. Rev. Genet. 7, 745. Excoffier, L., Laval, G., Schneider, S., 2005. Arlequin (version 3.0): an integrated software package for population genetics data analysis. Evol. Bioinform. Online. 1, 47. Excoffier, L., Lischer, H.E., 2010. Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Mol. Ecol. Resour. 10, 564–567. Felsenstein, J., 1985. Confidence limits on phylogenies: an approach using the bootstrap. Evolution 39, 783–791. Fielding, A.H., Bell, J.F., 1997. A review of methods for the assessment of prediction errors in conservation presence/absence models. Environ. Conserv. 24, 38–49. Fu, S.H., Fu, K.T., 1984. Crassulaceae. In: Chen, W.Q., Ruan, Y.Z. ( Eds.), Flora Reipublicae Popularis Sinicae, vol. 34. Science Press, Beijing, pp. 33–216. Fu, K.T., Ohba, H., 2001. Crassulaceae. In: Wu, C.Y., Raven, P.H. ( Eds.), Flora of China, vol. 8. Science Press and St. Louis: Missouri Botanical Garden Press, Beijing, pp. 202–268.

29

Gao, Q.B., Zhang, D.J., Duan, Y.Z., Zhang, F.Q., Li, Y.H., Fu, P.C., Chen, S.L., 2012. Intraspecific divergences of Rhodiola alsia (Crassulaceae) based on plastid DNA and internal transcribed spacer fragments. Bot. J. Linn. Soc. 168, 204–215. Guest, H.J., Allen, G.A., Carine, M., 2014. Geographical origins of North American Rhodiola (Crassulaceae) and phylogeography of the western roseroot, Rhodiola integrifolia. J. Biogeogr. 41, 1070–1080. Guindon, S., Gascuel, O., 2003. A simple, fast, and accurate algorithm to estimate large phylogenies by maximum likelihood. Syst. Biol. 52, 696–704. Guo, Z.T., Ruddiman, W.F., Hao, Q.Z., Wu, H.B., Qiao, Y.S., Zhu, R.X., Peng, S.Z., Wei, J.J. ,Yuan, B.Y., Liu, T.S., 2002. Onset of Asian desertification by 22 Myr ago inferred from loess deposits in China. Nature 416, 159–163. Hamilton, M.B., 1999. Four primer pairs for the amplification of chloroplast intergenic regions with intraspecific variation. Mol. Ecol. 8, 521–523. Harpending, H.C., 1994. Signature of ancient population growth in a low-resolution mitochondrial DNA mismatch distribution. Hum. Biol. 66, 591–600. Harris, N., 2006. The elevation history of the Tibetan Plateau and its implications for the Asian monsoon. Palaeogeogr. Palaeoclimatol. Palaeoecol. 241, 4–15. Hasumi, H., Emori, S., 2004. K-1 coupled GCM (MIROC) description. Center for Climate System Research, Univ. of Tokyo, Tokyo, Japan. Hewitt, G., 2000. The genetic legacy of the Quaternary ice ages. Nature 405, 907–913. 30

Hewitt, G.M., 2004. Genetic consequences of climatic oscillations in the Quaternary. Phil. Trans. R. Soc. B. 359, 183–195. Hijmans, R.J., Cameron, S.E., Parra, J.L., Jones, P.G., Jarvis, A., 2005. Very high resolution interpolated climate surfaces for global land areas. Int. J. Climatol. 25, 1965–1978. Hijmans, R.J., Guarino, L., Cruz, M., Rojas, E., 2001. Computer tools for spatial analysis of plant genetic resources data: 1. DIVA-GIS. Plant Genet. Resour. Newsl. 127, 15–19. Hughes, C.E., Atchison, G.W., 2015. The ubiquity of alpine plant radiations: from the Andes to the Hengduan Mountains. New Phytol. 207, 275–282. Li, J.J., Fang, X.M., Pan, B.T., Zhao, Z.J., Song, Y.G., 2001. Late cenozoic intensive uplift of Qinghai-Xizang plateau and its impacts on environments in surrounding area. Quat. Sci. 21, 381–391. Li, L., Abbott, R.J., Liu, B., Sun, Y., Li, L., Zou, J., Wang. X., Miehe, G., Liu, J., 2013. Pliocene intraspecific divergence and Plio-Pleistocene range expansions within Picea likiangensis (Lijiang spruce), a dominant forest tree of the Qinghai‐Tibet Plateau. Mol. Ecol. 22, 5237–5255. Librado, P., Rozas, J., 2009. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics 25, 1451–1452. Liu, L., Hao, Z.Z., Liu, Y.Y., Wei, X.X., Cun, Y.Z., Wang, X.Q., 2014. Phylogeography of Pinus armandii and its relatives: heterogeneous 31

contributions of geography and climate changes to the genetic differentiation and diversification of Chinese white pines. PLoS One 9, e85920. Luo, D., Yue, J.P., Sun, W.G., Xu, B., Li, Z.M., Comes, H.P., Sun, H., 2016. Evolutionary history of the subnival flora of the Himalaya‐Hengduan Mountains: first insights from comparative phylogeography of four perennial herbs. J. Biogeogr. 43, 31–43. Masters, B.C., Fan, V., Ross, H.A., 2011. Species delimitation - a geneious plugin for the exploration of species boundaries. Mol. Ecol. Res. 11, 154–157. Meng, L.H., Yang, R., Abbott, R.J., Miehe, G., Hu, T.H., Liu, J.Q., 2007. Mitochondrial and chloroplast phylogeography of Picea crassifolia Kom. (Pinaceae) in the Qinghai-Tibetan Plateau and adjacent highlands. Mol. Ecol. 16, 4128–4137. Mittermeier, R.A., Gil, P.R., Hoffman, M., Pilgrim, J., Brooks, T., Mittermeier, C.G., Lamoreux, J., Fonseca, G.A., 2005. Hotspots revisited: Earth's Biologically Richest and Most Endangered Terrestrial Ecoregions. Conservation International and Agrupacion Sierra Madre Press, Monterrey. Mittermeier, R.A., Turner, W.R., Larsen, F.W., Brooks, T.M., Gascon, C., 2011. Global biodiversity conservation: the critical role of hotspots. In: Biodiversity Hotspots. Springer, Berlin Heidelberg, pp. 3–22. Müller, K., 2005. SeqState. Appl. Bioinformatics 4, 65–69.

32

Myers, N., Mittermeier, R.A., Mittermeier, C.G., Fonseca, G.A.B., Kent, J., 2000 Biodiversity hotspots for conservation priorities. Nature 403, 853–858. Nordal, I., 1987. Tabula rasa after all? Botanical evidence for ice-free refugia in Scandinavia reviewed. J. Biogeogr. 14, 377–388. Ohba, H., 1978. Generic and infrageneric classfication of the old world Sedoideae (Crassulaceae). J. Fac. Sci. Univ. Tokyo, Bot. 3, 139–138. Ohba, H., 1987. Biogeography of the genus Rhodiola. Acta. Phytotax. Geobot. 38, 211–223. Opgenoorth, L., Vendramin, G.G., Mao, K., Miehe, G., Miehe, S., Liu, J.L., Ziegenhagen, B., 2010. Tree endurance on the Tibetan Plateau marks the world's highest known tree line of the Last Glacial Maximum. New Phytol. 186, 780–780. Otto-Bliesner, B.L., Marshall, S.J., Overpeck, J.T., Miller, G.H., Hu, A., Members CLIP, 2006. Simulating arctic climate warmth and icefield retreat in the last interglaciation. Science 311, 1751–1753. Peterson, A.T., Monica, P., Jorge, S., 2008. Rethinking receiver operating characteristic analysis applications in ecological niche modeling. Ecol. Model. 213, 63–72. Phillips, S.J., Dudík, M., 2008. Modeling of species distributions with Maxent: new extensions and a comprehensive evaluation. Ecography 31, 161–175.

33

Qin, A.L., Wang, M.M., Cun, Y.Z., Yang, F.S., Wang, S.S., Ran, J.H., Wang, X.Q., 2013. Phylogeographic evidence for a link of species divergence of Ephedra in the Qinghai-Tibetan Plateau and adjacent regions to the Miocene Asian aridification. PloS One 8, e56243. Rambaut, A., 2009. FigTree, version 1.3. 1. Computer program distributed by the author. http://treebioedacuk/software/figtree/. Rambaut, A., Drummond, A. J., 2007. Tracer version 1.4. Computer program and documentation distributed by the author. http://beast.bio.ed.ac.uk/ Tracer. Ren, G., Mateo, R.G., Liu, J.Q., Suchan, T., Alvarez, N., Guisan, A., Conti, E., Salamin, N., 2017. Genetic consequences of Quaternary climatic oscillations in the Himalayas: Primula tibetica as a case study based on restriction site-associated DNA sequencing. New Phytol. 213, 1500–1512. Renner, S.S., 2016. Available data point to a 4‐km‐high Tibetan Plateau by 40 Ma, but 100 molecular‐clock papers have linked supposed recent uplift to young node ages. J. Biogeogr. 43, 1479–1487. Richardson, J.E., Pennington, R.T., Pennington, T.D., Hollingsworth, P.M., 2001. Rapid diversification of a species-rich genus of neotropical rain forest trees. Science 293, 2242–2245. Rieseberg, L.H., Soltis, D.E., 1991. Phylogenetic consequences of cytoplasmic gene flow in plants. Evol. Trends Plants. 5, 65–84.

34

Rogers, A.R., 1995. Genetic evidence for a Pleistocene population explosion. Evolution 49, 608–915. Rogers, A.R., Harpending, H., 1992. Population-growth makes waves in the distribution of pairwise genetic-differences. Mol. Biol. Evol. 9, 552–569. Ronquist, F., Huelsenbeck, J.P., 2003. MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics 19, 1572–1574. Simmons, M.P., Ochoterena, H., 2000. Gaps as characters in sequence-based phylogenetic analyses. Syst. Biol. 49, 369–381. Spicer, R.A, Harris, N.B.W., Widdowson, M., Herman, A.B., Guo, S.X., Valdes, P.J., Wolfe, J.A., Kelley, S.P., 2003. Constant elevation of southern Tibet over the past 15 million years. Nature 421, 622–624. Spooner, D.M., 2016. Species delimitations in plants: lessons learned from potato taxonomy by a practicing taxonomist. J. Syst. Evol. 54, 191–203. Sun, Y.S., Ikeda, H., Wang, Y.J., Liu, J.Q., 2010. Phylogeography of Potentilla fruticosa (Rosaceae) in the Qinghai-Tibetan Plateau revisited: a reappraisal and new insights. Plant Ecol. Divers. 3, 249–257. Swofford, D.L., 2002. PAUP*. Phylogenetic Analysis Using Parsimony (*and Other Methods). Version 4.0b10. Sunderland, Sinauer Associates Press, Massachusetts.

35

Taberlet, P., Gielly, L., Pautou, G., Bouvet, J., 1999. Universal primers for amplification of three non-coding regions of chloroplast DNA. Plant Mol. Biol. 17, 1105–1109. Tamura, K., Stecher, G., Peterson, D., Filipski, A., Kumar, S., 1994. MEGA6: molecular evolutionary genetics analysis version 6.0. Bioinformatics 30, 2725. Wang, H., La, Q., Sun, K., Lu, F., Wang, Y., Song, Z., Wu, Q., Chen, J., Zhang, W., 2010. Phylogeographic structure of Hippophae tibetana (Elaeagnaceae) highlights the highest microrefugia and the rapid uplift of the Qinghai-Tibetan Plateau. Mol. Ecol. 19, 2964–2979. Wang, L.Y., Abbott, R.J., Zheng, W., Chen, P., Wang, Y.J., Liu, J.Q., 2009a. History and evolution of alpine plants endemic to the Qinghai-Tibetan Plateau: Aconitum gymnandrum (Ranunculaceae). Mol. Ecol. 18, 709–721. Wang, L.Y., Ikeda, H., Liu, T.L., Wang, Y.J., Liu, J.Q., 2009b. Repeated range expansion and glacial endurance of Potentilla glabra (Rosaceae) in the Qinghai-Tibetan Plateau. J. Integr. Plant Biol. 51, 698–706. Wen, J., Ickert-Bond, S.M., Appelhans, M.S., Dorr, L.J., Funk, V.A., 2015. Collections-based systematics: opportunities and outlook for 2050. J. Syst. Evol. 53, 477–488. Wen, J., Zhang, J.Q., Nie, Z.L., Zhong, Y., Sun, H., 2014. Evolutionary diversifications of plants on the Qinghai-Tibetan Plateau. Front. Genet. 5, 4.

36

White, T.J., Bruns, T.D., Lee, S.B., Taylor, J.W., 1990. Amplification and direct sequencing of fungal ribosomal RNA genes for phylogenetics. In: Innis, M., Gelfand, D., Sninsky, J., White, T. (Eds.), PCR Protocols: a Guide to Methods and Applications. Academic Press, San Diego, pp. 315–322. Wolfe, K.H., Li, W.H., Sharp, P.M., 1987. Rates of nucleotide substitution vary greatly among plant mitochondrial, chloroplast, and nuclear Dnas. Proc. Natl. Acad. Sci. U. S. A. 84, 9054–9058. Yamane, K., Yano, K., Kawahara, T., 2006. Pattern and rate of indel evolution inferred from whole chloroplast intergenic regions in sugarcane, maize and rice. DNA Res. 13, 197–204. Yang, F.S., Li, Y.F., Ding, X., Wang, X.Q., 2008. Extensive population expansion of Pedicularis longiflora (Orobanchaceae) on the Qinghai-Tibetan Plateau and its correlation with the Quaternary climate change. Mol. Ecol. 17, 5135–5145. Yi, T.S., Jin, G.H., Wen, J., 2015. Chloroplast capture and intra- and inter-continental biogeographic diversification in the Asian-New World disjunct plant genus Osmorhiza (Apiaceae). Mol. Phylogenet. Evol. 85, 10–21. Zhang, J.Q., Meng, S.Y., Allen, G.A., Wen, J., Rao, G.Y., 2014a. Rapid radiation and dispersal out of the Qinghai-Tibetan Plateau of an alpine plant lineage Rhodiola (Crassulaceae). Mol. Phylogenet. Evol. 77, 147–158.

37

Zhang, J.Q., Meng, S.Y., Wen, J., Rao, G.Y., 2014b. Phylogenetic relationships and character evolution of Rhodiola (Crassulaceae) based on nuclear ribosomal ITS and plastid trnL-F and psbA-trnH Sequences. Syst. Bot. 39, 441–451. Zhang, Q., Yang, R., Wang, Q., Liu, J.Q., 2005. Phylogeography of Juniperus przewalskii (Cupressaceae) inferred from the chloroplast DNA trnT-trnF sequence variation. Acta Phytotaxon. Sin. 43, 503–512. Zhang, Y.H., Volis, S., Sun, H., 2010. Chloroplast phylogeny and phylogeography of Stellera chamaejasme on the Qinghai-Tibet Plateau and in adjacent regions. Mol. Phylogenet. Evol. 57, 1162–1172. Zheng, B.X., Xu, Q.Q., Shen, Y.P., 2002. The relationship between climate change and Quaternary glacial cycles on the Qinghai-Tibetan Plateau: review and speculation. Quatern. Int. 97, 93–101. Zhou, W., Jin, J., Wu, J., Chen, H., Yang, J., Murphy, R.W., Che, J., 2017. Mountains too high and valleys too deep drive population structuring and demographics in a Qinghai-Tibetan Plateau frog Nanorana pleskei (Dicroglossidae). Ecol. Evol. 7, 240–252. Zuo, Y.J., Wen, J., Ma, J.S., Zhou, S.L., 2015. Evolutionary radiation of the Panax bipinnatifidus species complex (Araliaceae) in the Sino-Himalayan region of eastern Asia as inferred from AFLP analysis. Plant Syst. Evol. 53, 210–220.

38

Table 1 Genetic diversity and genetic differentiation of Rhodiola sect. Trifida at the species level and group level. HS

HT

GST

NST

All populations

0.314 (0.0425)

0.936 (0.0245)

0.665 (0.0443)

0.853 (0.0387)**

Clade I

0.539 (0.1803)

0.881 (0.0414)

0.388 (0.2237)

0.496 (0.2675)**

Clade II

0.178 (0.0840)

0.956 (0.0482)

0.814 (0.9140)

0.929 (0.0381)**

All populations

0.419 (0.0487)

0.959 (0.0077)

0.564 (0.5070)

0.637 (0.0587)ns

Clade A

0.364 (0.0527)

0.944 (0.0107)

0.714 (0.0555)

0.746 (0.0558)ns

Clade B

0.128 (0.0733)

0.765 (0.1150)

0.832 (0.0897)

0.897 (0.0684)ns

Clade D

0.055 (0.0545)

0.456 (0.1758)

0.980 (0.1184)

0.987 (.01740)ns

Datasets cpDNA

ITS

**NST is significantly different from GST (p < 0.01); ns, not significant. Table 2 Analysis of molecular variance (AMOVA) of chlorotypes and ITS ribotypes for Rhodiola sect. Trifida population cpDNA Source of variation

df

SS

ITS VC

PV (%)

F statistic 39

df

SS

VC

PV (%)

F statistic

Among groups

6

506.89

0.69928Va

19.30

FSC = 0.82419**

6

730.511

0.80637 Va

8.63

FSC = 0.66661**

Among populations

36

1152.68

2.40986 Vb

66.51

FST = 0.85812**

36

2390.968

5.69423 Vb

60.91

FST = 0.69537**

Within populations

547

272.45

0.51406 Vc

14.19

FCT = 0.01930**

450

1275.836

2.84785 Vc

30.48

FCT = 0.08626*

589

1932.31

3.62321

1.00

492

4397.316

9.34845

1.00

Total

df, degrees of freedom; SS, sum of squares; VC, variance components; PV, percentage of variation; FSC, correlation within populations relative to group; FST, correlation within populations relative to total; FCT, correlation within groups relative to total; ** P < 0.01, * P < 0.05, 10000 permutations.

40

Table 3 Results of the mismatch distribution analysis and neutrality tests of the two multiple-haplotype cpDNA clades and all populations Haplotype group

Tau

t (Mya)

SSD

p

Raggedness index

p

Tajima's D

p

FS

p

Clade I

8.503

-

0.421

0.105

0.115

0.063

1.7275

0.969

8.875

0.983

Clade II

6.301

-

0.558

0.062

0.164

0.018

1.896

0.971

1.477

0.759

All populations

0.743

0.135 (0.068-0.205)

0.005

0.699

0.106

0.854

-0.848

0.212

-12.709

0.032

41

Figure legends Fig. 1. A map showing the sites of sampled populations and the geographic distribution of haplotypes in sect. Trifida based on the cpDNA dataset. (A) Pie charts showing the proportion of haplotypes within each population. The numbers in the circles represented population number as in Table S1. Different shape on map represented different species, respectively. Dashed line on the map indicating the distribution area of sect. Trifida. (B) The locations of cpDNA private haplotypes throughout the distributional range. Each column referred to one population. Each box in a column referred to a single private haplotype. Fig. 2. A map showing the sites of sampled populations and the geographic distribution of ribotypes of sect. Trifida based on the ITS dataset. (A) Pie charts showing the proportion of ribotypes within each population. The numbers in the circles represented population number as in Table S1. Different shape on map represented different species, respectively. Dashed line on the map indicating the distribution area of sect. Trifida. (B) The locations of ITS private ribotypes throughout the distributional range. Each column referred to one population. Each box in a column referred to a single private ribotype. Fig. 3. Phylogenetic relationships based on cpDNA haplotypes. (A) The Bayesian tree topology of the 55 cpDNA haplotypes detected in sect. Trifida. Numbers above the branches are MP bootstrap support values (left) and Bayesian posterior probability (right). Different brunch color represents the three haplotype clades. (B) 42

NETWORK-derived genealogic relationships. The sizes of the circles in the network were proportional to the observed frequencies of the haplotypes in each population. Different colors representing different species. The red diamonds represented missing chlorotypes. Fig. 4. Phylogenetic relationships based on ITS ribotypes. (A) The Bayesian tree of the 47 ITS ribotypes detected in sect. Trifida. Numbers above the branches are MP bootstrap support values (left) and Bayesian posterior probability (right). (B) NETWORK-derived genealogic relationships. The sizes of the circles in the network are proportional to the observed frequencies of the haplotypes in each population. Different colors representing different species. The red diamonds representing missing ribotypes. Fig. 5. Divergence times of sect. Trifida based on ITS ribotypes estimated with BEAST. Blue bars indicating 95% highest posterior density intervals; numbers on the branch representing posterior probability and the number below representing mean age of the node. Fig. 6. Predicted distribution of sect. Trifida based on ecological niche modeling using Maxent during (A) the LIG. (B) the LGM-CCSM model. (C) the LGM-MIROC model. (D) present. White circles indicating localities that used for our simulation of the present study.

43

Supporting Information Fig. S1. Divergence times of sect. Trifida based on the plastid DNA haplotypes estimated with BEAST. Gray bars indicating 95% highest posterior density intervals. numbers on the branch representing posterior probability and the number below representing mean age of the node. Fig. S2. Historical demography for overall populations and in each regional group based on the plastid DNA dataset. Clade I and clade II corresponded to the Bayesian phylogenetic tree in Fig. 3A. Mismatch distribution showing histogram of observed mismatch frequencies and best-fit curve of the spatial expansion model. Table S1. Locations of populations of Rhodiola sect. Trifida sampled, sample sizes (N), frequencies of cpDNA haplotypes and ITS sequences per population, and estimates of haplotype diversity and nucleotide diversity for haplotypes and ribotypes within populations. Table S2. Haplotype composition of 43 sampled populations of Rhodiola sect. Trifida based on cpDNA dataset. Table S3. Ribotype composition of 43 sampled populations of Rhodiola sect. Trifida based on ITS dataset. Table S4. Bioclimatic variables (named BIO1 to BIO19). Variables marked with an asterisk (*) were used for the climatic niche models of Rhodiola sect. Trifida. Table S5. Information of the 203 localities used in the climatic niche models of Rhodiola sect. Trifida. 44

45

46

47

48

Graphical abstract

49

Highlights 1. The origin of sect. Trifida were in accordance with the last phase of the QTP uplifts 2. Multiple microrefugia existed on the QTP platform for sect. Trifida 3. The ecological niche modeling results also support the multiple microrefugia hypothesis

50