Rapid Intraspecific Diversification of the Alpine Species Saxifraga sinomontana (Saxifragaceae) in the Qinghai-Tibetan Plateau and Himalayas

An increasing number of phylogeographic studies have been conducted for plant species in the Qinghai-Tibetan Plateau (QTP) and its flanking mountains. However, these studies have mainly focused on the determination of glacial refugia and routes of inter-/post-glacial expansions. Rapid intraspecific diversification of plants in this region have not been thoroughly discussed. Herein, we investigate the effects of the Quaternary climate changes on population genetic structure and diversifications of a herbaceous alpine species, Saxifraga sinomontana, which may have an evolutionary time scale <5 million years in the QTP and Himalayan regions. Using a total of 350 individuals from 29 populations, we studied the evolutionary history of S. sinomontana by analyzing cpDNA trnL-trnF, rpl16 and nrDNA ITS sequences. A total of 89 haplotypes and 158 genotypes were detected for cpDNA and ITS sequences, respectively. Only a few haplotypes/genotypes were widespread, while an extremely large number of haplotypes/genotypes were restricted to single populations, which were scattered throughout the current geographical range of S. sinomontana. This suggests the existence of microrefugia of this species during the Quaternary glaciations. In addition, the relationships of the haplotypes/genotypes were almost completely not resolved by phylogenetic reconstruction. Combining characteristics in terms of high haplotype richness, large proportion of private haplotypes, and shallow haplotype divergence, we speculate that recent intraspecific diversification has occurred in S. sinomontana. Molecular clock analysis estimated that the onset diversification within S. sinomontana to be 1.09 Ma (95% HPD = 0.80–1.45), coinciding with the extensive Quaternary glaciations on the QTP which started ca. 1.17 Ma. The Quaternary climatic oscillations may have triggered rapid intraspecific diversification in this QTP-Himalayan species. However, large niche breadth, as well as introgression/hybridization between the studied species and its closely related sympatric saxifrages, may also played a role to some extent on the current genetic structure of S. sinomontana, which need to be further studied.

An increasing number of phylogeographic studies have been conducted for plant species in the Qinghai-Tibetan Plateau (QTP) and its flanking mountains. However, these studies have mainly focused on the determination of glacial refugia and routes of inter-/post-glacial expansions. Rapid intraspecific diversification of plants in this region have not been thoroughly discussed. Herein, we investigate the effects of the Quaternary climate changes on population genetic structure and diversifications of a herbaceous alpine species, Saxifraga sinomontana, which may have an evolutionary time scale <5 million years in the QTP and Himalayan regions. Using a total of 350 individuals from 29 populations, we studied the evolutionary history of S. sinomontana by analyzing cpDNA trnL-trnF, rpl16 and nrDNA ITS sequences. A total of 89 haplotypes and 158 genotypes were detected for cpDNA and ITS sequences, respectively. Only a few haplotypes/genotypes were widespread, while an extremely large number of haplotypes/genotypes were restricted to single populations, which were scattered throughout the current geographical range of S. sinomontana. This suggests the existence of microrefugia of this species during the Quaternary glaciations. In addition, the relationships of the haplotypes/genotypes were almost completely not resolved by phylogenetic reconstruction. Combining characteristics in terms of high haplotype richness, large proportion of private haplotypes, and shallow haplotype divergence, we speculate that recent intraspecific diversification has occurred in S. sinomontana. Molecular clock analysis estimated that the onset diversification within S. sinomontana to be 1.09 Ma (95% HPD = 0.80-1.45), coinciding with the extensive Quaternary glaciations on the QTP which started ca. 1.17 Ma. The Quaternary climatic oscillations may have triggered rapid intraspecific diversification in this QTP-Himalayan species. However, large niche breadth, as well as introgression/hybridization between the studied species and its closely related sympatric saxifrages, may also played a role to some extent on the current genetic structure of S. sinomontana, which need to be further studied.

INTRODUCTION
Evolutionary diversifications are usually considered as one of the main mechanisms that accumulate high levels of plant biodiversity in mountainous regions (Hughes and Eastwood, 2006;Pennington et al., 2010;Hughes and Atchison, 2015). This is the case for alpine plants in the Qinghai-Tibetan Plateau (QTP) and its surrounding Himalayas and Hengduan Mountains region (HHM), from where a large number of rapid diversifications have been recorded (see Wen et al., 2014 and references therein). The driving factors that have triggered such diversifications could be extrinsic (e.g., orogeny and climate change) (Hoorn et al., 2010), intrinsic (e.g., innovation of novel traits, polyploidization, hybridization, and niche shifts), or a combination (Ebersbach et al., 2017b). Among these factors, orogenic events are proposed to have played a disproportionate contribution to plant diversifications in the QTP and HHM at generic, even higher taxonomic levels (e.g., Liu et al., 2006;Zhang et al., 2012Zhang et al., , 2014Favre et al., 2015Favre et al., , 2016Ebersbach et al., 2017a,b;Xing and Ree, 2017). In contrast to deeply divergent lineages in which diversifications were mainly triggered by ancient orogeny, intraspecific differentiations of recently divergent species could be shaped by extrinsic events at a small timescale, such as the Quaternary glaciations, which started ca. 1.17 million years ago (Ma; Zheng et al., 2002). Repeated environment alternations between arid/cold glacial and humid/warm inter-glacial episodes may have caused habitat fragmentation which promoted isolation of populations, and thereby facilitated intraspecific differentiation (Hickerson et al., 2010;Gao et al., 2012Gao et al., , 2016Jia et al., 2012). Due to extremely complex topolography , extensive width of ecological niches (Ni and Herzschuh, 2011), lack of unified glaciers (Zheng et al., 2002), as well as effects of monsoon systems (Zhang et al., 2000), the evolutionary histories of plants in the QTP and HHM during the Quaternary glaciations seem to be more complicated than in North America and Europe (see Qiu et al., 2011 and references therein). Species that have different distribution ranges, population sizes, mating systems, life forms, and bio-characteristics may have experienced distinctive patterns of divergence and evolutionary history during glaciations (Shahzad et al., 2017). The accumulation of phylogeographic data for as many species as possible can provide a more complete picture of the intraspecific evolutionary history of plants in the QTP and HHM.
Encompassing ca. 500 species in at least 13 sections (Pan et al., 2001;Gao et al., 2015Gao et al., , 2017Tkach et al., 2015;Ebersbach et al., 2017b), Saxifraga L. is the largest genus of Saxifragaceae s.str. (Soltis, 2007). Its distribution is primarily concentrated in the mountainous and arctic regions of the Northern Hemisphere. This species-rich genus has widely been employed in the fields of systematics and phylogeography to reveal patterns and processes of plant diversification in arctic and alpine regions (e.g., Conti et al., 1999;Vargas, 2000Vargas, , 2001Abbott and Comes, 2003;Oliver et al., 2006;Westergaard et al., 2010;DeChaine et al., 2013;Gao et al., 2015Gao et al., , 2017Ebersbach et al., 2017a,b). The QTP and HHM is a biodiversity center of Saxifraga, harboring nearly half the species, mainly represented by the two species-rich sections Ciliatae Haworth and Porphyrion Tausch (Pan et al., 2001;Ebersbach et al., 2017a). Although immigration, recent rapid radiation, as well as lineage persistence should all contribute to the extremely high diversity and endemism of Saxifraga in the QTP and HHM (Ebersbach et al., 2017a), rapid geographic and adaptive radiations have proven to be disproportionally important for the QTP-HHM lineages of sects. Ciliatae and Porphyrion (Ebersbach et al., 2017b). However, nearly all of these studies have focused on the roles of diversifications or radiations on species richness of Saxifraga in the QTP and HHM, while case studies of intraspecific differentiations are still scarce for this species-rich Saxifraga region.
Saxifraga sinomontana J. T. Pan & Gornall is one of the common Saxifraga species in the QTP and HHM, inhabiting shrublands, alpine/marshy meadows and rock crevices at elevations of between 2,700 and 5,300 m a.s.l. (Pan et al., 2001). This perennial herb is extraordinarily variable in morphology, consecutively from small individuals (<5 cm) with a single flower to tall ones (more than 20 cm) with multiple flowers. This high variation probably suggests potential intraspecific diversification in this species. Phylogenetic analyses have supported a position of S. sinomontana in sect. Ciliatae subsect. Hirculoideae Engl. & Irmsch. (Gao et al., 2015;Tkach et al., 2015). It is characterized by pedicels with brown crisped villi as well as erect sepals with crisped villi at margins and abaxial surface (Pan et al., 2001). However, species relationships within clade of S. subsect. Hirculoideae were not well resolved in previous studies (Gao et al., 2015;Tkach et al., 2015), partly due to recent rapid radiations triggered by uplifts of the Hengduan Mountains (Ebersbach et al., 2017a,b). In fact, it seems that nearly all of the ca. 110 species in S. subsect. Hirculoideae diverged <5 Ma (Ebersbach et al., 2017a,b), or even more recently (ca. 2 Ma; Gao et al., 2015). As a member of this subsection, intraspecific differentiation of S. sinomontana must have occurred more recently, and probably affected by climatic oscillations during the Quaternary glaciations. Therefore, it is suitable for revealing associations between intraspecific diversification and the Quaternary climatic oscillations.
In the present study, we sequenced two chloroplast DNA (cpDNA) fragments (trnL-trnF and rpl16) and nuclear ribosomal DNA internal transcribed spacer (ITS) to reveal the population genetic structure and evolutionary history of S. sinomontana. In particular, we want to address: (1) whether intraspecific diversification has occurred in this QTP-HHM species; and (2) whether climatic oscillations during the Quaternary glaciations promoted the intraspecific diversification if it occurred.

Sample Collection
Fresh leaves of one to 24 individuals (according to the population size) were sampled from each population, spaced at least 5 m apart. In total, leaf materials of 350 individuals from 29 populations were collected across the extant distribution range of S. sinomontana (Table 1, Figure 1). Leaves were dried in silica gel. Voucher specimens of all populations are deposited in the

DNA Extraction, PCR Amplification, and Sequencing
Total genomic DNA was extracted from silica-dried leaves using the modified CTAB method (Doyle and Doyle, 1987). For polymerase chain reactions (PCR), primers of "c" and "f " for trnL-trnF (Taberlet et al., 1991), "F71" and "R1516" for rpl16 (Kelchner and Clark, 1997), and "18S" and "26A" for ITS (Wen and Zimmer, 1996) were employed for amplification. PCR reactions were carried out in a total volume of 50 µl containing 0.4 µl (1.5 units) of Taq polymerase, 5.0 µl of 10 × PCR buffer (with Mg 2+ ), 2.0 µl of 10 mM dNTPs, 1.0 µl of 10 pM of each primer and 2.0 µl (10-20 ng) total genomic DNA. The PCR profile included an initial pretreatment of 10 min at 94 • C, followed by 30 cycles of 1 min denaturation at 94 • C, 50 s at primer-specific annealing temperatures, 1 min elongation at 72 • C, and a final extension at 72 • C for 10 min. Annealing temperatures were 49 • C for trnL-trnF, 58 • C for rpl16, and 55 • C for ITS. PCR products were checked on 1% agarose gels and then purified using a CASpure PCR Purification Kit (CASarray, Shanghai, China). Purified PCR products were sequenced in both directions with the primers used for amplification on an ABI PRISM 3730xl genetic analyzer.

Haplotype Isolation
The chromatograms of each sequence were contrasted for accuracy via visual inspection with Chromas ver. 2.6.4 (available at http://www.technelysium.com.au). Primer regions were removed and DNA sequences were aligned in MEGA ver. 7.0.26 (Kumar et al., 2016) with minor subsequent adjustments. The ragged tails of the alignments were trimmed to ensure an uniform ending. All cpDNA sequences were identified to different haplotypes using DnaSp ver. 5.10 (Librado and Rozas, 2009). The two cpDNA fragments (trnL-trnF and rpl16) were then concatenated into a single matrix for subsequent analyses. For ITS sequences, a site was identified as heterozygous when a double peak occurred in the same position of both strands, with the weakest signal reaching at least a quarter of the strength of the strongest (Fuertes Aguilar et al., 1999;Fuertes Aguilar and Nieto Feliner, 2003). The haplotype phase for ITS sequences was reconstructed using PHASE 2.1 (Stephens et al., 2001;Stephens and Donelly, 2003) as implemented in DnaSP. Subsequent analyses of the ITS dataset were based on isolated haplotypic sequence types. All newly generated sequences in this study have been deposited in GenBank (accession no. MH432703-MH433033).

Phylogenetic Analysis
Phylogenetic relationships of the detected cpDNA/ITS haplotypes were reconstructed by means of Maximum Parsimony (MP), Maximum Likelihood (ML), and Bayesian Inference (BI). The MP analyses were conducted in PAUP 4.10b (Swofford, 2002) with all characters and transitions/transversions equally weighted. A heuristic search with 100 random-taxonaddition replicates was performed with ACCTRAN character optimization, MULPARS + TBR branch swapping, MulTrees and STEEPEST DESCENT options selected. Estimates of bootstrap support (BS) were calculated using 1,000 replicates with the above settings. jModelTest ver. 2.1.4 (Darriba et al., 2012) was used to choose the best-fit substitution models for the ML and BI analyses according to the Akaike Information Criteria (AIC). Model of GTR+I+G was selected for the cpDNA and ITS datasets. The ML analyses were conducted using RAxML ver. 8.1.21 (Stamatakis, 2014) as implemented in raxmlGUI ver. 1.5b2 (Silvestro and Michalak, 2012) with a selection of ML + rapid bootstrap and support assessment using 1,000 rapid bootstraps. The BI analyses were performed using MrBayes ver. 3.2.6 (Ronquist and Huelsenbeck, 2003;Ronquist et al., 2012) and the same substitution modes as the ML analyses. Two simultaneous Markov Chain Monte Carlo (MCMC) analyses were run for 10 million generations, sampling every 1,000 generations with the first 25% trees discarded as burn-in. Due to limited parsimony informative sites among cpDNA/ITS haplotypes, phylogenetic reconstructions of  Table 1. Names of haplotypes and genotypes are the same as in Table 2. MP, ML, and BI resulted in an almost complete lack of resolution in haplotype relationships (Supplementary Files 1, 2). We have noticed that high proportion of variable sites were represented as indels in the aligned cpDNA dataset. These indels may contain important phylogenetic information, and their involvement, to some extent, may improve phylogenetic resolution of cpDNA dataset. Indels were then coded using method described by Simmons and Ochoterena (2000) as implemented in FastGap ver. 1.2 (available at: http://www.aubot.dk/FastGap_home.htm). Phylogenetic reconstructions of indels-coded cpDNA haplotypes were repeated as described above, however, phylogenetic resolution was not yet improved (Supplementary File 1). To detect genealogical relationships among sequences with shallow genetic divergences, cpDNA/ITS haplotype networks were constructed by the program NETWORK ver. 5.0.0.3 (available at: http://www.fluxus-engineering.com), using median-joining method and MP calculation (Polzin and Daneshmand, 2003).

Population Genetic Structure Analysis
Average gene diversity within populations (H S ), total gene diversity (H T ), as well as two coefficients of population genetic differentiation, G ST and N ST , were estimated as described by Pons and Petit (1996) using the program PERMUT (available at: http:// www.pierroton.inra.fr/genetics/labo/Software/PermutSSR). To assess whether the overall population differentiation exhibited a phylogeographic structure, a comparison was made between G ST and N ST using a permutation test with 1,000 permutations. Significantly larger N ST value compared with G ST value suggests the presence of phylogeographic structure (Pons and Petit, 1996). Populations with less than three individuals were discarded from the PERMUT analysis, namely populations DR (one individual), LZI (two individuals), and SD (two individuals). Analysis of molecular variance (AMOVA) was performed using Arlequin ver. 3.5.2 (Excoffier and Lischer, 2010) to partition genetic variation within and among populations. Gene diversity (h) and nucleotide diversity (π) of each population were calculated using Arlequin ver. 3.5.2 (Excoffier and Lischer, 2010). We further calculated haplotype/genotype diversity at the species level as an indicator of variation by dividing the number of haplotypes/genotypes recovered by the number of individuals assayed (Oliver et al., 2006). The number of private haplotypes (N ph ) and genotypes (N pg )  Population codes are the same as those given in Table 1. within each population was also calculated to assess genetic diversity.

Population Demographic Analysis
Mismatch distribution analysis based on the cpDNA dataset was conducted to infer the demographic history of S. sinomontana in the overall populations. The shape of the frequency graph of pairwise differences is expected to be multimodal in samples drawn from populations at demographic equilibrium, whereas it is usually unimodal in populations having passed through a recent demographic or range expansion with high levels of migration between neighboring demes (Slatkin and Hudson, 1991;Rogers and Harpending, 1992;Harpending et al., 1998;Ray et al., 2003;Excoffier, 2004). One thousand parametric bootstrap replicates were used to generate an expected distribution under a model of sudden demographic expansion. The sum of squared deviations (SSD) was used as a statistical test to accept or reject the hypothesis of sudden population expansion. Harpending's raggedness index (HRI, Harpending, 1994) and its significance were calculated to quantify the smoothness of the observed mismatch distribution. Neutrality tests of Tajima's D (Tajima, 1989) and Fu's F S (Fu, 1997) based on the cpDNA dataset were also carried out to assess possible expansion. Large negative values of Tajima's D and Fu's F S statistics should be observed under the expansion hypothesis. All of these demographic tests were implemented using Arlequin ver. 3.5.2 (Excoffier and Lischer, 2010).

Divergence Time Estimation
The ITS dataset was used to estimate divergence times of detected haplotypes using BEAST ver. 1.8.4 (Drummond et al., 2012). Since ancient fossil calibrations might bias divergence time estimation toward younger ages in tip clades (Sauquet et al., 2012;Gao et al., 2017), intraspecific divergence time of S. sinomontana was estimated using a substitution rate of 5.03 × 10 −9 substitutions per site per year (s/s/y) of ITS fragment as estimated by Ebersbach et al. (2017a) in Saxifragaceae. Nucleotide substitution model of GTR+I+G, a strict molecular clock and the Yule tree prior were applied. Two independent MCMC analyses were run for 10 million generations, with sampling every 1,000 generations. The maximum clade credibility tree was summarized in TreeAnnotator ver. 1.8.4 (Drummond et al., 2012) with 25% burn-in.

Haplotype/Genotype Distribution and Phylogenetic Relationships
The alignment length of sequences of individual cpDNA fragments was 897 bp and 821 bp for trnL-trnF and rpl16, respectively. However, there was an adenine homoploymer at site 308 and a tyrosine homopolymer at site 355 in trnL-trnF, to ensure that variable sites were actual mutations instead of sequencing errors, nucleotides of and between homopolymers were deleted in the following analyses of this intergenic spacer. In addition, adenine and tyrosine homopolymers were also detected at sites 125 and 722, respectively, in rpl16, we deleted nucleotides of the two homopolymers for this fragment in the following analyses. Based on the concatenated cpDNA sequences (trnL-trnF and rpl16), 89 haplotypes (H1-H89) were identified ( Table 2). Sequence lengths of these haplotypes varied from 1,303 to 1,420 bp, with an alignment length of 1,603 bp. Only a few haplotypes, e.g., H18 (occurred in 35 individuals from 12 populations), H17 (31 individuals from eight populations), H15 (22 individuals from nine populations), H32 (17 individuals from six populations), and H10 (10 individuals from six populations), were widespread across the extant distribution range of S. sinomontana. However, large amounts of haplotypes (63 out of 89) were private, i.e., restricted to single populations. Most of these exhibited extremely low level of occurrence frequency ( Table 2).
With respect to ITS dataset, a total of 158 ITS genotypes were discovered from the 350 individuals of S. sinomontana ( Table 2), and their sequence length varied from 717 to 718 bp with an alignment length of 718 bp. Distribution pattern of the detected ITS genotypes was similar to those of cpDNA haplotypes, that is, a small number of widespread genotypes (e.g., T1/T6, T6/T6) vs. a high proportion (88%) of population-private genotypes. Based on the 158 different ITS genotypes, we isolated 116 haplotypes (T1-T116) using the program PHASE 2.1. The following analyses for ITS dataset were based on the isolated haplotypes.
Since some private haplotypes/genotypes were identified as singletons, i.e., represented by single individuals, which can be putatively caused by amplification and/or sequencing errors, reduplication of amplification and sequencing was carried out for such individuals to omit the extremely low error rates of methodology.
Phylogenetic reconstructions for both cpDNA and ITS haplotypes detected in S. sinomontana, by means of MP, ML, and BI, showed a "comb-like" topology. This resulted in an almost complete lack of resolution of haplotype relationships (Supplementary File 1). Median-joining networks were then conducted for these shallow-divergent cpDNA/ITS haplotypes. In both cpDNA and ITS networks, haplotypes with high distribution frequency (e.g., H10, H17, and H18 for cpDNA dataset; T1, T3, and T6 for ITS) were located in the central positions of individual networks, while population-specific haplotyes with low frequency generally occupied network tips (Figures 2, 3). In addition, the ITS haplotype network showed a star-shaped phylogeny (Figure 3). Mutation steps among neighboring cpDNA haplotypes ranged from one to nine, with the highest divergence occurring between H87 and H11/H22 by nine mutation steps (Figure 2). However, divergence among neighboring ITS haplotypes was even shallower, generally discriminated by one mutation step (Figure 3).

Population Genetic Diversity and Phylogeographic Structure
Gene diversity (h) of cpDNA dataset within the 29 populations ranged from 0 to 1.0000, and nucleotide diversity (π) from 0 to 0.000425, with mean values of 0.8064 and 0.000144, respectively ( Table 2). Gene diversity and nucleotide diversity based on ITS data among the 29 populations ranged from 0.2048 to 1.0000, and 0.000590 to 0.008123, with average values of 0.7577 and 0.003001, respectively ( Table 2). Haplotype and genotype diversity at species level were 0.25 and 0.45, respectively. Nearly all detected populations fixed private cpDNA haplotypes and/or ITS genotypes. Both cpDNA and ITS datasets revealed high values of total genetic diversity (H T = 0.978 for cpDNA, and 0.896 for ITS), as well as average withinpopulation diversity (H S = 0.784 for cpDNA, and 0.747 for ITS), suggesting a disproportionate contribution of within-population diversity to the total in S. sinomontana. This was confirmed by analysis of molecular variance (AMOVA), which showed that 77.94 and 60.80% of total genetic variation was found within populations based on cpDNA and ITS datasets, respectively (Table 3)

Population Demographic Analyses and Divergence Time
Both neutrality tests and mismatch distribution analysis based on cpDNA dataset suggested recent range or demographic expansion of S. sinomontana. Neutrality tests of Tajima' D (−2.134, P < 0.001) and Fu' F S (−4.923, P < 0.001) statistics produced significantly negative values for the overall gene pool, suggesting recent expansion across the distribution range of S. sinomontana. This was supported by mismatch distribution analysis, in which unimodal was drawn from the overall populations (Figure 4). Besides, both indices of sum of squared deviations (SSD, 0.004, P > 0.05) and Harpending's raggedness index (HRI, 0.004, P > 0.05) were not significant, showing no deviation of observed mismatch distribution from simulation under a model of sudden demographic expansion.
The ITS dataset was employed to estimate the onset time of intraspecific divergence in S. sinomontana. Based on a mean substitution rate of 5.03 × 10 −9 substitutions per site per year (s/s/y) of ITS fragments as estimated by Ebersbach et al. (2017a) in Saxifragaceae, the onset of intraspecific divergence of S. sinomontana was estimated to be 1.09 Ma (95% HPD = 0.80-1.45), corresponding to the Quaternary climatic oscillations (Figure 5).

DISCUSSION
An increasing number of plant phylogeographic studies have been conducted in the QTP and its flanking mountains (e.g., Zhang et al., 2005;Meng et al., 2007;Opgenoorth et al., 2010;Wang et al., 2010; also see review by Qiu et al., 2011). However, the aims of these studies mainly focused on the determination of glacial refugia and routes of inter-/postglacial expansions. Rapid intraspecific differentiations of plants in this region have not been well discussed, although rapid evolutionary diversifications have been reported in an increasing number of QTP plant groups (e.g., Liu et al., 2006;Zhang et al., 2012Zhang et al., , 2014Favre et al., 2015Favre et al., , 2016    and Ree, 2017; also see review by Wen et al., 2014), including the largest circum-QTP section Ciliatae of Saxifraga (Ebersbach et al., 2017b), in which S. sinomontana was embedded. Genetic structure of extant populations of S. sinomontana differs from that of infra-generic circumarctic-alpine species S. oppositifolia L. (Abbott and Comes, 2003) and S. hirculus L. (Oliver et al., 2006). However, genetic structure of S. sinomontana is comparable to that of sympatric S. tangutica (Gengji et al., 2018).

Rapid Intraspecific Diversification in S. sinomontana
Recently, high geological heterogeneity and rapid climate oscillations have been highlighted as potential forces for species diversification in alpine plants (Simões et al., 2016;Ebersbach et al., 2017b). Rapid plant diversification could result in phylogenetic uncertainty (Richardson et al., 2001;Hughes and Eastwood, 2006;Wen et al., 2013;Ebersbach et al., 2017b) which is attributable to shallow genetic divergence among species because of limited accumulation of mutations. In addition, lineages that have experienced recent and rapid diversifications usually exhibit high level of regional endemism (Liu et al., 2006;Zhang et al., 2012Zhang et al., , 2014Favre et al., 2015Favre et al., , 2016, probably due to short migration history, geographic barriers and a lack of suitable but unoccupied habitats. Similarly, if rapid intraspecific diversifications have occurred in a plant species, large amounts but shallow divergence of haplotypes, as well as a high proportion of private haplotypes should be expected, especially for species with limited gene flow among populations. In this study, an extremely large number of cpDNA haplotypes (89) and ITS genotypes (158) are identified across the 350 individuals from 29 populations of S. sinomontana. The haplotype and genotype diversity (the number of haplotypes/genotypes divided by the number of individuals) at the species level of S. sinomontana are 0.25 and 0.45, respectively. These are higher than the values for other QTP herbaceous species, e.g., Aconitum gymnandrum (cpDNA haplotype diversity/ITS genotype diversity: 0.04/not available; , Pedicularis longiflora (0.03/not available; Yang et al., 2008), Rhodiola alsia (0.20/0.08; Gao et al., 2012), and Meconopsis integrifolia (0.03/0.17; Yang et al., 2012). In fact, S. sinomontana seems to show the highest haplotype/genotype diversity among the studied plant species in the QTP and adjacent regions to date, although the comparison is somewhat arbitrary since different DNA fragments have been employed for different species. A high proportion of these detected haplotypes/genotypes can be described as private (71% of cpDNA haplotypes and 88% of ITS genotypes), i.e., observed in single populations. Most of these haplotypes/genotypes exhibited extremely low frequencies, with some even represented as singletons. However, limited mutations have been detected among such a large number of haplotypes, and most neighboring haplotypes are separated by limited mutational steps. This is even more apparent in the ITS haplotype network, with most differing by only a single site mutation. Limited informative mutation sites among both cpDNA and ITS haplotypes resulted in an almost complete lack of resolution of haplotype relationships as reconstructed by means of MP, ML, and BI (Supplementary Files 1, 2), indicating recent rapid intraspecific differentiation in S. sinomontana. Shallow divergence of detected cpDNA haplotypes can be further confirmed by low level of mean nucleotide diversity (0.000144) across assayed populations of S. sinomontana. In addition, it should be notable that ITS dataset revealed a relatively high level of mean nucleotide diversity (0.003948), which could be attributed to high amounts of hybridization sites. Considering haplotype richness, proportion of private haplotypes, and haplotype divergence, it is reasonable to speculate that recent intraspecific diversification occurred in this QTP-Himalayan species, S. sinomontana. Global cooling and climatic oscillations, which started in the Middle Miocene and intensified continuously until the end of the Pleistocene (Miao et al., 2012;Favre et al., 2015), are considered as having played a crucial role for QTP plant diversifications Xing and Ree, 2017;Ebersbach et al., 2018). Intraspecific diversification of S. sinomontana may also have been triggered by climatic changes but should be more recent, e.g., glacial and interglacial intervals in the Quaternary. Our molecular clock analysis estimated that the onset of diversification within S. sinomontana to be 1.09 Ma (95% HPD = 0.80-1.45), coinciding with the first of the four extensive Quaternary glaciations on the QTP which started ca. 1.17 Ma (Zheng et al., 2002). It should kept in mind that this time estimation was based on the mean ITS substitution rate of the large family Saxifragaceae (Ebersbach et al., 2017a), and it is probably an under-estimate since the diversification rate in clade Ciliatae subsect. Hirculoideae, which includes S. sinomontana, was detected to be ∼2.5 times higher than the background rate (Ebersbach et al., 2017b). This means that the onset differentiation time of S. sinomontana may be more recent, falling well within the drastic climate oscillations of the Quaternary. It is likely that repeated developments and retreats of glaciers during Quaternary fragmented distribution range of S. sinomontana into isolated patches, finally facilitating rapid intraspecific radiation.
However, other factors may have also played a role to some extent in the context of rapid intraspecific diversification of S. sinomontana. Among them, attention should be paid to large niche breadth, which was considered to be one of the two main factors that drove rapid diversification in sect. Ciliatae subsect. Hirculoideae (Ebersbach et al., 2018). S. sinomontana is an extraordinarily variable species (Pan et al., 2001). Populations that inhabit humid environments at relatively low elevations, such as shrublands and alpine/marshy meadows, usually exhibit tall flowering stems (>20 cm) and bear many flowers. While populations that inhabit rock crevices and calcareous rocks near nival line exhibit loose tufts with short flowering stems (<5 cm) and solitary flower. Intermediates and reticulate variations occur between the two extreme phenotypes. The extraordinary variation in morphology of S. sinomontana appears to be environmentally induced (Oliver et al., 2006), and it probably indicates a large niche breadth of this species. Species with a large niche breadth usually possess wide environmental tolerance, and this enables a broad distribution range (Warren et al., 2010), which could increase the opportunity for diversification by: (1) increased chance of population separation and differentiation (Gehrke and Linder, 2011;Tanentzap et al., 2015;Matuszak et al., 2016); (2) increased ecological opportunity (Losos, 2010); and (3) low chance of extinction during climatic changes (Ebersbach et al., 2018).
Overall, the differentiation among populations of S. sinomontana is low as revealed by both cpDNA and ITS datasets (G ST = 0.199 and 0.166, F ST = 0.2206 and 0.3920 for cpDNA and ITS, respectively). One possible explanation should be efficient gene flow among populations (Brochmann et al., 2003;Oliver et al., 2006). However, a large proportion of private haplotypes/genotypes have been detected in S. sinomontana, and almost all populations harbor private haplotypes/genotypes, indicating limited gene flow among populations. Besides, the seeds of S. sinomontana have no special adaptations for dispersal. Although mechanisms of pollination and seed dispersal are not clear in S. sinomontana, they are probably comparable to those of S. hirculus, a species showing similar morphology and close phylogenetic relationships (Pan et al., 2001;Gao et al., 2015). Seeds of S. hirculus are mainly dispersed by dropping near the parent plants with an average dispersal distance of 13 cm (Olesen and Warncke, 1989). Furthermore, although S. hirculus retains partial characteristics of outcrossing, self-pollination seems to be more important for reproductive assurance in rough weather conditions because of lack of effective pollinators (Li et al., 2013). Thus, effective gene flow is not a plausible explanation for the low level of population differentiation in S. sinomontana. However, this pattern could be caused by recent rapid intraspecific diversification of S. sinomontana, which is characterized by large numbers of private haplotypes/genotypes vs. few widespread ones. High frequencies of a small number of widespread haplotypes/genotypes across populations of S. sinomontana could decrease population differentiation. However, a high proportion of private haplotypes/genotypes could increase within-population diversity but contribute little to among-population differentiation. This pattern could be produced by their extremely low frequency, scattered distribution and shallow divergence. This is also explicable for the detection of lack of phylogeographic structure across the distribution range of S. sinomontana, as also have been observed in Potentilla glabra , Stellera chamaejasme (Zhang et al., 2010), Rhodiola alsia (Gao et al., 2012), and Rhodiola chrysanthemifolia (Gao et al., 2016).
A word of caution needs to be injected here because there may be other players in this story. In particular, we do not yet know to what extent, if any, hybridization and/or introgression play a role on extant genetic structure of S. sinomontana. According to our field investigation, S. sinomontana can be locally sympatric to its close relatives S. heleonastes H. Smith, S. congestiflora Engl. & Irmsch., S. tangutica L., and S. pseudohirculus Engl., even rarely sympatric to S. przewalskii Engl. and S. tibetica Losinsk. Sometimes, individuals of S. sinomontana can be <10 cm apart from species mentioned above, which provides a geographic opportunity for hybridization and/or introgression. Thus, considering the revealed genetic structure of S. sinomontana and rapid radiation of S. sect. Ciliatae subsect. Hirculoideae (Gao et al., 2015;Ebersbach et al., 2017b), hybridization and/or introgression between S. sinomontana and its sympatric relatives cannot be ruled out. A combining study of local S. sinomontana populations and its close relatives in sympatry can give us a complete picture of evolution of this alpine plant complex.

Microrefugia for S. sinomontana During Glaciations
Phylogeographic patterns of plant species on the QTP during the Quaternary glaciations have been reviewed (Qiu et al., 2011), and their corresponding genetic structures as reflected in current populations have also been throughly discussed (Gao et al., 2016). In this study, the decline in genetic diversity from southeastern edge of the QTP to the plateau platform was not detected (e.g., Zhang et al., 2005;Meng et al., 2007), and regional differentiation centers could not be identified (e.g., Gao et al., 2012). Instead, large numbers of private haplotypes/genotypes are scattered across the distribution range of S. sinomontana, and populations with high genetic diversity show an even distribution, suggesting existence of microrefugia of S. sinomontana on the QTP during the Quaternary glaciations. In fact, considering high heterogeneity of topography of the QTP, especially of the Hengduan Mountains , as well as the absence of unified glaciers (Zheng et al., 2002), it is of great opportunity to provide suitable micro-environments for cold-tolerant herbs (Gao et al., 2016), even shrubs and trees (Opgenoorth et al., 2010;Wang et al., 2010), to survive glaciations in situ. However, a recent expansion signal was detected across the distribution range of S. sinomontana by both neutrality tests and mismatch distribution analysis based on cpDNA dataset. As we discussed above, the current genetic structure of S. sinomontana, in conjunction with the fact of low dispersal ability of seeds, provides evidence against extensive horizontal range expansion across the distribution range of S. sinomontana. Thus, the detected expansion signal probably represents demographic expansion or altitudinal migration in response to repeated glacier developments and retreats, as suggested in Rhodiola chrysanthemifolia (Gao et al., 2016), Potentilla glabra , and P. fruticose (Shimono et al., 2010). It is likely that population of S. sinomontana was continuous in the QTP-HHM region before the Quaternary glaciations, with some widespread haplotypes/genotypes (e.g., H15, H17, H18 for cpDNA; T1/T6, T6/T6 for ITS) across its distribution range. The following developments and retreats of glaciers during the Quaternary may have fragmented distribution range of S. sinomontana into isolated patches, finally facilitating in situ allopatric divergence. The ancient genetic structure of S. sinomontana may have erased to some extent due to bottleneck effects and genetic drifts during the Quaternary glaciations, which was then replaced by current genetic pattern represented by large numbers of private haplotypes/genotypes.