Quaternary climate change and habitat preference shaped the genetic differentiation and phylogeography of Rhodiola sect. Prainia in the southern Qinghai–Tibetan Plateau

Abstract There are two long‐standing biogeographic hypotheses regarding the glacial survival of plant species in the Qinghai–Tibetan Plateau (QTP): the in situ survival hypothesis and the tabula rasa hypothesis. We tested these two hypotheses in a phylogeographic study of Rhodiola sect. Prainia, a monophyletic section with ecologically divergent lineages. Molecular data from the nuclear internal transcribed spacer, six plastid markers and 13 nuclear microsatellite loci were analyzed for 240 individuals from 19 populations of this section. Environmental data were used to analyze the niches of major phylogenetic lineages within this section and to model changes in their distributions since the Last Glacial Maximum (LGM). We found that Rhodiola sect. Prainia consists of three evolutionary lineages: all populations of R. stapfii, R. prainii populations at the southern edge of the QTP, and R. prainii populations in the interior part of the QTP. During the LGM, the survival of R. prainii in the interior part of the QTP corresponded with the in situ survival hypothesis, while R. stapfii most probably survived the LGM in a manner corresponding with the tabula rasa hypothesis. The evolutionary history of different lineages of this section was shaped by topography, climate change, and lineage‐specific habitat preferences.

patterns of alpine species, that is, the two abovementioned hypotheses, influence of topography and climate changes on the distribution of alpine species should be considered, with a particular focus on taxon-specific characteristics (Papadopoulou & Knowles, 2016).
The Qinghai-Tibetan Plateau (QTP) consists of many mountains, flatlands, and valleys, and it provides a wide range of potential refugia and geographical barriers for testing hypotheses regarding glacial survival. The timing and extent of glaciations in the QTP during the Last Glacial Maximum (LGM) were still disputable, but it is generally agreed that there was no ice sheet covering the entire plateau and the valleys of the Yarlung Zangbo River were mostly unglaciated (Owen & Dortch, 2014). Fossil pollen records in the QTP reveal that vegetation zones have moved and expanded from the southeast to the northwest, from the edge of the plateau to the hinterland, and from low elevations to high elevations during the Holocene (Hou, Yang, Cao, Chongyi, & Wang, 2015). These changes in plant distribution allowed both postglacial colonization from peripheral areas and in situ survival via vertical movement along mountain slopes. Many alpine trees and shrubs occurring in the QTP, including Potentilla fruticosa (Li, Shimono, Shen, & Tang, 2010), Potentilla glabra (Wang, Ikeda, Liu, Wang & Liu 2009), the Juniperus tibetica complex (Opgenoorth et al., 2010), Spiraea alpina (Zhang et al., 2012), and a lineage of Hippophae tibetana , probably survived in situ during Quaternary glaciations. In contrast, some other alpine trees and shrubs, including Juniperus przewalskii (Li, Zhang, Liu, Källman, & Lascoux, 2011;Zhang, Chiang, George, Liu, & Abbott, 2005), Picea crassifolia (Meng et al., 2007), and Tsuga dumosa (Cun & Wang, 2010), migrated through these barriers and recolonized the QTP. The herb species of the QTP show similar differences in their demographic histories. For example, four subnival herbs survived locally in the southern QTP (Luo et al., 2016), while the alpine herbs Metagentiana striata (Chen et al., 2008) and Pedicularis longiflora (Yang, Li, Ding, & Wang, 2008) recolonized the QTP from its southeastern edge.
Previous studies of various taxa in the QTP provide data supporting both of the hypotheses mentioned above, but there seems to be no fixed pattern regarding the evolutionary trajectories followed by plant species in the QTP. Additionally, no species has been reported to recolonize the QTP from its southern edge (Qu, Lei, Zhang, & Lu, 2010;Yu et al., 2017), which is unexpected given the reports of southern refugia in Europe and North America (Hewitt, 2000), the presence of several south-north valleys connecting the southern edge of the QTP to its interior (Garzione, DeCelles, Hodkinson, Ojha, & Upreti, 2003), and the distributions of certain species that occur on both sides of the Himalayas (Opgenoorth et al., 2010;Qu et al., 2010;Ren et al., 2017).
The southern QTP is a rugged region characterized by high biodiversity and high endemism (Zhang, Ye, & Sun, 2016), and its most prominent topographic feature is the Yarlung Zangbo River valleys between the Himalayas and the Gangdese-Nyainqentanglha Range.
Steep slopes within 10 km of the Yarlung Zangbo River and its tributaries lead to differences in elevation of more than 3,000 m (Wang et al., 2014) and significant altitudinal variation in vegetation, which ranges from dry riparian vegetation up to alpine meadows (Wu & Wu, 1998). Opgenoorth et al. (2010) argued that the topography of the southern QTP and the neighboring Hengduan Mountains is more varied in comparison with that of the northern QTP, and their diverse habitats enabled more species to survive glaciations in situ. In contrast, several species have recolonized the southern QTP from refugia in the southeastern edge of the QTP (Cun & Wang, 2010;Yang et al., 2008;Yu et al., 2017). Thus, differences in the demographic histories of species in the southern QTP may have been caused by their specific biological characteristics such as habitat preferences or dispersal capacity.
To test whether different habitat preferences can lead to different evolutionary trajectories under a shared climatic and geological background, and whether plant species can recolonize the QTP from its southern edge, we focused on Rhodiola sect. Prainia, a monophyletic taxon consisting of Rhodiola stapfii, and R. prainii. Currently, these two species have a similar distribution pattern in the southern QTP, but they occur in different habitats.
R. stapfii and R. prainii are both distributed disjunctively on mountains around the valleys of the Yarlung Zangbo River and the southern edge of the QTP. Moreover, they form a species pair (Zhang, Meng, Allen, Wen, & Rao, 2014), and share many features.
The most conspicuous difference between these species is their distinct habitat preferences. R. prainii occurs on rock slopes or cliffs, while R. stapfii is restricted to alpine meadows. Alpine meadows are common at high elevations in the QTP, while rock slopes are distributed discontinuously along river valleys or suture zones. If the biogeographic histories of R. stapfii and R. prainii were dominated by shared geographic events and climate change, then a synchronic evolutionary history would be expected. Otherwise, the influence of different habitat preferences on glacial survival and postglacial dispersal should be taken into consideration.
In this study, 19 populations of R. stapfii and R. prainii were sampled across their distribution areas, including both the southern edge and interior of the QTP. Plastid sequences, nuclear internal transcribed spacer (ITS) sequences, and nuclear microsatellites (simple sequence repeats, SSRs) were used to elucidate the phylogenetic relationships, phylogeographic patterns, and demographic histories of the studied species. In addition, their historical distributions and migration routes were reconstructed using species distribution models (SDMs). Furthermore, we assessed whether these two closely related species behaved differently under the shared influence of Quaternary climate change, and whether they recolonized the interior of the QTP from its southern edge.

| Sampling and DNA extraction
Samples were collected from the mountains and shallow valleys of the southern edge and interior of the QTP (Figure 1). In this study, the southern edge of the QTP refers to the southward-flowing watersheds at the southern slope of the Himalayas. Leaves were collected from 240 individuals from 19 populations of R. prainii and R. stapfii (eight R. stapfii populations and eleven R. prainii populations) in the field throughout their distributions and dried in silica gel (see Table S1). Although a detailed field investigation was conducted at the southern edge of the QTP, only one population was sampled in this region for R. prainii (Pra11) and R. stapfii Samples were also collected from closely related species as outgroups (Table S1)

| Genetic marker design and detection
Microsatellite markers and informative plastid DNA markers were developed from high-throughput sequencing data as follows (Table S2).
Five plastid genomes (three R. prainii and two R. stapfii) were acquired following Dong, Xu, Cheng, Lin, and Zhou (2013). Five pairs of plastid primers were designed to cover the most informative regions on the plastid genomes according to the following criteria: (a) the PCR product should be less than 500 bp in length; (b) primers should be located in coding regions; (c) the fragment should cover as much variation as possible; (d) SSR or poly A/T regions should be avoided (Table   S4). With regard to traditional plastid markers, only trnH-psbA (Sang, Crawford, & Stuessy, 1997) was used because it contains abundant variation in R. stapfii and R. stapfii (Table S3). Thus, six plastid markers in total were used in this study. As to nuclear markers, nuclear reads were assembled de novo into short contigs, from which SSR loci were detected and their primers designed. Next, 165 pairs of SSR primers were chosen randomly and screened for suitable stability and polymorphisms in a sample consisting of three individuals from each population. Finally, 13 robust and highly informative SSR loci were selected for genotyping of all individuals of both species (Table S5). Traditional marker ITS sequences were also used (Mayuzumi & Ohba, 2004). PCR was performed following Zhang et al. (2014). SSRs were amplified and detected following Schuelke (2000).

| Phylogenetic analysis
To confirm the relatedness of the two studied species, a phylogenetic analysis including several closely related species was carried out. DNA sequences were aligned with muscle 3.8.31 (Edgar, 2004) and checked manually. ITS ribotypes and concatenated plastid haplotypes were generated with dNAsp 5.10.01 (Librado & Rozas, 2009). Haplotype networks of concatenated plastid sequences were drawn using the TCS method implemented in popArt 1.7 (Leigh & Bryant, 2015). Phylogenetic analyses were conducted within a Bayesian framework using beAst 1.8.4 (Drummond & Rambaut, 2007). The ITS tree was dated by setting the divergence time between Rhodiola and Phedimus in beAst to 21.02 Ma with a standard deviation of 6.5 Ma . Under the Bayesian information criterion (Schwarz, 1978), jmodeltest 2.1.7 (Darriba, Taboada, Doallo, & Posada, 2012) selected K2 + G (Kimura, 1980) as the best-fit substitution model for the ITS sequences, whereas TN93 (Tamura & Nei, 1993) was selected for the concatenated plastid sequences. Convergence was checked using TRACER 1.7 (Rambaut, Drummond, Xie, Baele, & Suchard, 2018). Phylogenetic analyses were also conducted using the maximum parsimony (MP) method implemented in pAup 4 (Swofford, 2003) and the Bayesian approach in mrbAyes 3.2.6 (Ronquist et al., 2012; Figure S1-S4).

| Population genetic structure
For SSR data, diversity indices (H e , H o ) were calculated in ArlequiN 3.5.2.2 (Excoffier & Lischer, 2010). The genetic structure was inferred based on SSR loci using structure 2.3.4 (Pritchard, Stephens, & Donnelly, 2000) assuming an admixture model. The burn-in was set to 20,000, and 50,000 generations were run after it. The K-value was then set to range from 1 to 12, and the computation was repeated 10 times for each Kvalue. Convergence was reached if alpha value arrived at a plateau in each run and 10 different runs with the same K-value resulted in similar results (Porras-Hurtado et al., 2013). The ΔK criterion implemented in structure hArvester v0.6.94 (Earl & vonHoldt, 2012) could not provide a clear best K-value in this study ( Figure S5; Table S5), possibly because the ΔK criterion is insensitive to the second-level genetic structure under complex schemes (Evanno, Regnaut, & Goudet, 2005). A clear population structure emerged when K = 3 or 4. clumpp 1.1.2 (Jakobsson & Rosenberg, 2007) and distruct 1.1 (Rosenberg, 2004) were used to summarize and visualize the structure results.
In addition, the genetic differentiation among populations within each evolutionary lineage was estimated using analysis of molecular variance (AMOVA) in ArlequiN.

| Recent population dynamics
Neutrality tests on DNA sequences and bottleNeck analysis on SSR data were used to test for possible postglacial expansions of the two studied species. dNAsp was used to calculate Tajima's D, Fu and Li's D*, and Fu and Li's F* based on ITS sequences to identify signatures of demographic expansions (Librado & Rozas, 2009), and 10,000 coalescent simulations were run to estimate the p-value of each index. Changes in population size for different groups were calculated using a two-phase model (TPM) with 95% single-step mutations and a variance among multiple steps of 12 in bottleNeck 1.2 using the SSR data (Piry, Luikart, & Cornuet, 1999). A sign test, a standardized differences test, and Wilcoxon's signed rank test were performed in bottleNeck to determine whether there was a heterozygosity deficit. Frequency of private SSR alleles (alleles detected in only one population of Sect. Prainia; Slatkin & Takahata, 1985) in each population were calculated in R (for scripts see Appendix S2 files). Ten individuals were randomly selected 100 times for the target population to correct the uneven population size, and the number of private alleles was averaged across replicates. The approximate Bayesian computation (ABC) approach implemented in diyAbc 2.1.0 (Cornuet et al., 2014) was used to compare five possible evolutionary scenarios for populations of R. stapfii and R. prainii ( Figure 3). All populations were divided into four groups before ABC according to the structure results at k = 4 ( Figure 4). Mean number of alleles and mean genic diversity of each group are summarized in diyAbc. Finally, 5,000,000 simulated data sets were generated, and the 10,000 data sets that were closest to the observed data set were used to estimate each parameter.

| Environmental factor analysis and niche modeling
To understand how the distributions of R. prainii and R. stapfii changed after the LGM, they were modeled under three climatic scenarios from the CCSM4 model (current climate at 30 arc-seconds resolution, mid-Holocene 6,000 years ago and the LGM 21,000 years ago at 2.5 arc-minutes resolution) with mAxeNt 3.3.3k (Phillips & Dudík, 2008).
Nineteen ecological factors were acquired from Worldclim 1.4 (available at: www.world clim.org; Hijmans, Cameron, Parra, Jones, & Jarvis, 2005). Current environmental data for 19 sampling sites as well as 30 occurrence records from our field investigation and online data sets (Table S7) were used in this study. To avoid collinearity, Pearson pairwise correlation analysis of environmental factors was conducted, and one factor was eliminated in each pair with a correlation value higher than 0.8 (Table S8; Ren et al., 2017). The selected factors were also used to describe the niche of each population with principal component analysis (PCA) using the "princomp" function in r 3.3.1 (R Core Team, 2016). SDMs were built under the maximum entropy method implemented in mAxeNt with the settings reported by Papeş and Gaubert (2007). Records of the occurrence of R. prainii at the southern edge of the QTP were not included in the SDM analysis because they form a lineage different from other R. prainii populations in the interior QTP ( Figure 5) and occur in different environments. Selected rasters of environmental factors were cropped to span from 26°N to 33°N and from 78°E to 100°E before modeling.

| Microsatellite marker design and selection
According to analyses of high-throughput sequencing data from five individuals of the Rhodiola sect. Prainia species (two R. stapfii and three R. prainii), an average of 15,499,029 clean reads were acquired for each individual, from which 2,514 pairs of nuclear SSR primers were designed. Finally, 13 robust and highly informative SSR loci were used for the population genetic analysis of the Sect. Prainia species: R. stapfii and R. prainii (Table S5).

| Phylogenetic and haplotype analysis
The final set of DNA sequences consisted of 570 bp of the ITS and 2081 bp of concatenated plastid sequences. Fourteen ITS ribotypes and nine plastid haplotypes were detected across R. stapfii and R. prainii (Table 1). The phylogenetic trees for the ITS and concatenated plastid sequences were incongruent (Figures 5, 6). The ITS tree showed that all populations of Rhodiola sect. Prainia con-   (Table 2), which is consistent with the low level of differentiation found within each cluster by the structure analysis at K = 3 (Figure 4). S t a 1 S t a 2 S t a 3 S t a 4 S t a 5 S t a 6 S t a 7 S t a 8 P r a 1 P r a 2 P r a 3 P r a 4 P r a 5 P r a 6 P r a 7 P r a 8 P r a 9 P r a 1 0 P r a 1 1

| Recent population dynamics
were no contradictions because "recent" in bottleNeck was defined as within approximately the past 2N e -4N e generations (Piry et al., 1999), while the results of the pooled population analyses reflected their common demographic history, which was much longer than that of a single population (Table 1). ABC analysis identified scenario 3 as the best evolutionary model for Rhodiola sect. Prainia, which suggests that populations of R. stapfii diverged first, after which populations of R. prainii split into two lineages on both sides of the Himalayas (Figure 3). Considering that the generation time of the studied species is estimated to be approximately five years based on our field observation and previous studies (Galambosi, 2014;Wu, Shang, Dai, & Yan, 2001), the split time between populations Pra1 to Pra2 and Pra3 to Pra10 (T 1 ) was estimated at ~511 [305-729] generations (Table 4), which should be after the LGM.
The split time between the interior plateau (Pra1 through Pra10) and the southern edge populations (Pra11) of R. prainii (T 2 ) was estimated at ~9,320 [3,980-27,400] generations, (Table 4), which would be earlier than the LGM.

| Environmental factor analysis and SDMs
To reduce multicollinearity, one factor of each pair of environmental factors was eliminated when the correlation coefficient of the pair was >0.8. Thus, only four factors were left: annual mean temperature, isothermality, annual precipitation, and precipitation seasonality (Table S8). PCA of the four environmental factors at all known locations of Rhodiola sect. Prainia showed that the first two principal components accounted for 79.6% of the total variance (Table 5) (17) Hap5 (2) Hap6 (4) Hap7 (3) Chap5 (12) Chap6 (1)  eastern and western parts of the southern QTP and the area along the middle Himalayas (Figure 8c).

| Phylogenetic relationship among lineages of Rhodiola sect. Prainia
The incongruence of the topologies of the gene trees constructed from the ITS and plastid sequences may have been the result of inadequate phylogenetic information, inaccurate species limits, incomplete lineage sorting, unrecognized paralogy, or hybridization (Funk & Omland, 2003). The high posterior probabilities of the three major clades on the Rhodiola sect. Prainia ITS tree ar-

| Evolutionary histories of R. prainii and R. stapfii
The lack of genetic structure within R. stapfii suggests that current populations could be derived from only one refugium after the LGM (Table 2;  hypothesis. The Himalayas served as westward migration corridors rather than geographic barrier for R. stapfii, which is a recurrent pattern exhibited by other alpine species of this region (Wallis, Waters, Upton, & Craw, 2016;Yu et al., 2017). In addition, the valleys of the Yarlung Zangbo River were also potential refugia for R. stapfii (Figure 8a), in which case the tabula rasa hypothesis would be inaccurate. Evolutionary scenarios in mountainous regions could be far more complex than the simple dichotomy of the in situ survival versus tabula rasa hypotheses (Sersic et al., 2011;Stehlik, 2003), and a mixed survival scenario could not be ruled out in the case of R. stapfii.
ABC analysis showed that the establishment of the three lineages of Rhodiola sect. Prainia was earlier than the LGM (Table 4; Figure 3), and this finding was also supported by the dated phylogenetic tree of the Rhodiola sect. Prainia ITS sequences ( Figure 5). These results suggest that the three lineages survived several earlier glaciations that were more extensive than the LGM (Lehmkuhl & Owen, 2005).
Therefore, due to its limited distribution and strong contractions during the LGM, the evolutionary history of Sect. Prainia before the LGM cannot be conjectured based on the results of this study.

| Valley refugia
Both R. prainii and R. stapfii survived glaciations in deep river valleys. During glaciations, the valleys of the Yarlung Zangbo River were mostly unglaciated (Owen & Dortch, 2014) and provided glacial refugia for many local species at a wide range of elevations because of their significant climatic buffering capacity (Frenzel, Bräuning, & Adamczyk, 2003;Liang, He, Jia, Sun, & Chen, 2017;Opgenoorth et al., 2010;Zhang, Comes, & Sun, 2011). The significant role of valleys as glacial refugia is also documented in Patagonia (Sersic et al., 2011), indicating that "valley refugia" might be a common pattern in mountainous regions.
However, deep valleys are not sufficient as glacial refugia during the LGM. The current distribution of R. stapfii covers several deep river valleys (Figure 1), but the species had only one refugium during the LGM. The SDM results also suggest that southward valleys along the southern edge of the QTP were unsuitable for R. stapfii and R.

| Lack of recolonization of the interior from the southern edge of the QTP
In some cases, southward valleys at the southern edge of the QTP could serve as glacial refugia, but they do not act as sources TA B L E 2 Analysis of molecular variance (AMOVA) for the SSR data of interglacial recolonization (Ren et al., 2017;Opgenoorth et al., 2010;Yu et al., 2017). Under the dominating influence of the Indian monsoon, eastward or southeastward Himalayan valleys showed significant differences in precipitation patterns along elevation gradients in comparison with those of southward valleys (Bookhagen & Burbank, 2010). As a consequence, lineages occurring at the southern edge of the QTP have diverged genetically and ecologically from interior populations for millions of years, which has left them incapable of adapting to the current environments in the interior part of the QTP, as observed for population Pra11 of R. prainii (Figures 2, 7).
Another possible explanation could be that recolonization from the southern edge of the QTP was impeded by populations that survived in situ in the interior QTP. Genetic discrepancies resulting from isolated glacial refugia could be maintained even under postglacial contact and introgression (Hewitt, 1996;Wallis et al., 2016).
If a large number of populations survived in situ in the interior QTP, migration from the southern edge would lead to a genetic signature of introgression instead of a genetic pattern of recolonization from the southern edge. Such introgression was not found in Rhodiola sect. Prainia, but it has been documented in the Juniperus tibetica complex (Opgenoorth et al., 2010), suggesting that preexisting populations in the interior QTP could block recolonization from the southern edge in some cases. This situation also exists in the southeastern QTP, where some cold-adapted species recolonized the Himalayas from the Hengduan Mountains across biogeographic barriers, but other cold-adapted species survived in situ (Cun & Wang, 2010;Li et al., , 2011Liu et al., 2013;Opgenoorth et al., 2010;Yu et al., 2015).

| Migration barriers and corridors
The Previous studies on the demographic history of species in the southern QTP identified two common westward migration routes: (a) along the Himalayas and (b) along the Yarlung Zangbo River valleys (Qiu et al., 2011;Yu et al., 2017). In this study, R. stapfii migrated along the Himalayas, whereas R. prainii expanded in both directions along the valleys of the Yarlung Zangbo River, indicating that the migratory pattern of a species is influenced by its habitat preference, as manifested by characteristics such as habitat preferences, in addition to the locations of glacial refugia.
Although the SDMs predicted that there were suitable habitats for R. stapfii in the southwestern Himalayas ( Figure 8c) and alpine meadows are widespread in this region (Wang et al., 2016)

| Population size dynamics of herbs in the QTP
Although many cold-adapted species have experienced glacial expansions and interglacial contractions (Stewart et al., 2010), phylogeographic studies of alpine herbs on the southern QTP have not consistently reached the same conclusions (Wang, Abbott, et al., 2009). Four subnival herbs, which occur mostly in rocky environments at elevations between 3,400 and 5,000 m, experienced downslope expansion during the LGM and upslope contraction during postglaciation periods (Luo et al., 2016). Another herb species, Pedicularis longiflora, which is distributed in wet meadows and along hill streams at elevations between 2,600 and 5,300 m, experienced contraction during glacial advancement and expansion during glacial retreats (Yang et al., 2008). Asynchronous population size changes were documented in Primula tibetica, another herb with four genetic lineages occurring in wet meadows and along streams at elevations between 2,600 and 5,000 m. Two lineages of P. tibetica experienced a scenario of "'expansion-shrinkage-expansion'," while two other lineages experienced a scenario of "'expansion-shrinkage"' (Ren et al., 2017). In this study, R. stapfii and R. prainii of the inner QTP were sampled from elevations of 4,100-5,000 m. R. stapfii from alpine meadows experienced strong contraction during the LGM and expansion afterward, which is similar to the pattern displayed by Pedicularis longiflora (Yang et al., 2008). However, R. prainii of the inner QTP showed changes similar to those experienced by meadow herbs during glaciations instead of those experienced by species in rocky environments. One possible explanation for this finding could be that R. prainii of the inner QTP had a narrow adaptation amplitude that restricted the distribution of this lineage during glaciations.
In comparison with other alpine species from rocky environments, R. prainii is usually found on rocks with a relatively large amount of moss (per. obs.), which suggests that micro-habitat preferences might be another explanation for the migratory pattern of this species. These findings suggest that any explanations of differences in changes in population size for herbs on the southern QTP are inevitably influenced by taxon-specific adaptation characteristics besides cold tolerance.

| CON CLUS ION
This study showed that Rhodiola sect. Prainia, a taxon endemic to the southern QTP, consists of three evolutionary lineages rather than two species as previously described. R. prainii distributed at the southern edge of the QTP should be treated as an independent phylogenetic species (Nixon & Wheeler, 1990). The survival strategy of R. prainii in the interior part of the QTP during the LGM corresponded with the in situ survival hypothesis, while that of R. stapfii probably corresponded most closely with the tabula rasa hypothesis.
Taken together with previous studies, the findings of this study suggest that southward valleys at the southern edge of the QTP are either not appropriate glacial refugia, or not sources of recolonization to the interior QTP (Cun & Wang, 2010;Opgenoorth et al., 2010;Ren et al., 2017;Yang et al., 2008;Yu et al., 2015).

ACK N OWLED G M ENTS
This work is supported by the National Natural Science Foundation of China (Grant No. 31470313). We are grateful to Jian-Qiang Zhang for his help with the field work, Yun-Cheng Duan for laboratory work, and Wen-Pan Dong for data analysis.

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R S CO NTR I B UTI O N
G.R. and Z.W. conceived the ideas. Z.W. and S.M. conducted the field research. Z.W. did the laboratory work and analyzed data, supervised by G.R. All authors participated in drafting the manuscript.

DATA ACCE SS I B I LIT Y
DNA sequences: GenBank accessions MG917775-MG917968.