Mountains too high and valleys too deep drive population structuring and demographics in a Qinghai–Tibetan Plateau frog Nanorana pleskei (Dicroglossidae)

Abstract Pleistocene glacial–interglacial climatic oscillations greatly shaped the current genetic structure of many species. However, geographic features may influence the impact of climatic cycling. Distinct geographic and environmental characters between northern and southern parts of the eastern Qinghai–Tibetan Plateau (EQTP) facilitate explorations into the impacts of geographic features on species. The northern parts of EQTP contain large areas of marsh, and the environment is rather homogeneous. In contrast, the southern EQTP harbors complex alpine valleys and a much more heterogeneous setting. We evaluate DNA sequence variation from both the mitochondrial and nuclear genomes in Nanorana pleskei, a species endemic to the EQTP. Hypothesis testing on the evolutionary history of N. pleskei indicates that northern populations can disperse freely, but alpine valleys isolate southern populations. Demographic histories between northern and southern populations also differ. Northern populations appear to have experienced population expansions, while southern frogs exhibit a far more stable demographic history. By combining climatic analyses and species' distribution models, our study suggests that geographic and environmental features drive the differences between the northern and southern EQTP.

tions into the impacts of geographic features on species. The northern parts of EQTP contain large areas of marsh, and the environment is rather homogeneous. In contrast, the southern EQTP harbors complex alpine valleys and a much more heterogeneous setting. We evaluate DNA sequence variation from both the mitochondrial and nuclear genomes in Nanorana pleskei, a species endemic to the EQTP. Hypothesis testing on the evolutionary history of N. pleskei indicates that northern populations can disperse freely, but alpine valleys isolate southern populations. Demographic histories between northern and southern populations also differ. Northern populations appear to have experienced population expansions, while southern frogs exhibit a far more stable demographic history. By combining climatic analyses and species' distribution models, our study suggests that geographic and environmental features drive the differences between the northern and southern EQTP.

K E Y W O R D S
amphibians, demographic history, eastern Qinghai-Tibetan plateau, Pleistocene climatic oscillations, population structure

| INTRODUCTION
Pleistocene climatic cycling is one of the most important drivers of contemporary diversity in many temperate species and communities. Following the global cyclical cooling-warming events in the Quaternary, the distributions of various organisms experienced concomitant expansions and contractions. During this process, vicariance drove the divergence of populations, and the disappearance of barriers resulted in dispersal and secondary contact. These cyclical events appear to have shaped the current genetic structure of most temperate species (Hewitt, 2000(Hewitt, , 2004. Notwithstanding, geographic features also influence the distributions of species, especially in regions not covered by a unified ice sheet. Geographic barriers may promote vicariance in association with climatic changes. For example, barriers may stop dispersal and, thus, restrict expansion. Complex geographic features may harbor suitable microhabitats in the form of microrefugia. These refugia may buffer the impacts of climate changes and shape patterns of genetic diversity by stabilizing populations. In contrast, fluctuations in population size of organisms are easy to happen in the homogeneous environment. The complex geographic features and dramatic changes in climate make the eastern edge of the Qinghai-Tibetan Plateau (EQTP) an ideal place for exploring the impacts of the environment on organisms.
The EQTP consists of distinctive northern and southern parts based on their differing environments. The northern EQTP (nEQTP) occurs east of the Bayanhar Mountains and north of the Qionglai Mountains ( Figures S1 and 1). It largely consists of high-elevation plains, and the environment is more homogeneous than in the southern EQTP (sEQTP). The nEQTP has a few rolling hills void of steep slopes. The Yellow River flows through it (Figures 1 and S1). The upstream Dadu River contains two major branches-the Maerke and Duoke riversthat flow from north to south. Large areas of marsh form along these rivers.
The sEQTP exists south of the Bayanhar Mountains (You & Yang, 2013). This region belongs to the Hengduan Mountains, which are famous for their complex environmental features and diverse ecozones. These features contribute to it being a biodiversity hot spot (Myers, Mittermeier, Mittermeier, Da Fonseca, & Kent, 2000). The sEQTP consists of heterogeneous alpine valleys made by four longitudinal mountains separated by several major branches of the Yangtze River. From east to west, the Dadu River flows between the Qionglai and Daxueshan mountains, the Yalong River flows between the Daxueshan and Shaluli mountains, and the Jinsha River flows west of Shaluli Mountain ( Figures S1 and 1). These rivers form deeply incised valleys. The environment, including temperature, precipitation, and vegetation forms, changes drastically in the valley. The environment is arid and hot at bottom of these valleys, especially in the middle of Hengduan Mountains which range from 28°N to 32°N (You & Yang, 2013). The environments in the valleys potentially form barriers to dispersal. Drastic environmental differences between sEQTP and its adjacent region, nEQTP, provide opportunities to test how the environment influences genetic diversity and further understand the factors driving the regional biodiversity forming.
Although much research focuses on species in this region, very few studies evaluate the effects of rivers as geographic barriers on organisms. Among the exceptions, Scutiger boulengeri is an Alpine stream frog of the sEQTP, and its matrilineal genealogy suggests that the Yalong and Dadu rivers are effective barriers (Li, Chen, Tu, & Fu, 2009). In the lesser striped shrew (Sorex bedfordiae), divergences of the matrilines roughly correspond to large rivers (Chen et al., 2015).
In contrast, rivers alone failed to explain the population structure of Rana kukunoris (Zhao, Dai, & Fu, 2009) and a range-wide study of this species , which did not test specifically for river effects, also failed to discover that rivers served as barriers to dispersal.
The EQTP experienced drastic climatic changes during the Quaternary. Although no unified ice sheet formed, isolated montane glaciers, ice caps, and valley glaciers existed, especially in the sEQTP (Zhou, Wang, Wang, & Xu, 2006). Pollen records suggest that permafrost/desert steppe largely covered the QTP during the Last Glacial Maximum (LGM), yet forest islands existed in the southeastern and eastern ranges (Tang & Shen, 1996;Tang, Shen, Kong, Wang, & Liu, 1998). Glacial refugia, which have contributed greatly to the formation of genetic diversity pattern and strongly impacted demographic F I G U R E 1 Sampling sites and the maternal genealogy of Nanorana pleskei based on trees yield from beast. Site numbers detailed in Table  S1. Colors of lineage correspond with localities on the map. Bootstrap proportions (maximum-likelihood [ML] and maximum parsimony [MP] trees) ≥70% and Bayesian posterior probabilities ≥95% were treated as strong support ( ). Bootstrap proportions ≥50% and Bayesian posterior probabilities ≥75% were treated as weakly support ( ). Bootstrap proportions <50% and Bayesian posterior probabilities <75% were treated as not support ( ). Support rates of the nodes display following Bayesian inference/ML/MP. Numbers in nodes indicate divergence times yield from beast. The two dash-line boxes indicate nEQTP and sEQTP. The black lines indicate rivers   (Chen et al., 2008;Li et al., 2011;Wang et al., 2011;Yue & Sun, 2014). Refugia in the nEQTP are also reported (Jia, Liu, Wang, Zhou, & Liu, 2011;Zhou et al., 2013).
Because geological features may influence demographic history, and given the differences in geography between northern EQTP and southern EQTP, we test hypotheses related to the impact of climatic cycling and geography on organisms.
To reveal the impacts of climatic changes and geographic features on organisms, we use an endemic, medium-sized frog, Nanorana pleskei (Dicroglossidae). It occurs commonly in open habitats of the EQTP at elevations ranging from 3,000 to 4,500 m ( Figure S1). The species prefers a lentic environment, including marshes, pools, and ponds (Fei, Hu, Ye, & Huang, 2009

| Sampling and sequencing
We sampled 181 individuals from 21 localities covering most of the distribution ( Figure S1) of N. pleskei ( Figure 1; Table S1). Two samples of N. parkeri were included as the out-group taxa based on previous studies (Che et al., 2010). All collecting adhered to animal-use protocols approved by the Kunming Institute of Zoology Animal Care and Ethics Committee.
DNA was extracted from tissue samples using the standard phenol-chloroform extraction protocol (Sambrook & Russell, 2001).
Fragments were subjected to purification using EXO-SAP (Takara Bio Inc., Shiga, Japan). Both forward and reverse directions were sequenced to preclude sequencing errors.

| Genealogy reconstruction and population relationship inference
We used Bayesian inference (BI), maximum-likelihood (ML), and maximum parsimony (MP) to infer the matrilineal genealogy from mtDNA sequences. The best partition strategy and nucleotide substitution model for each data partition were determined using partItIOnfInder 1.1.1 (Lanfear, Calcott, Ho, & Guindon, 2012). After assigning the corresponding model to each partition, BI analyses were performed using Mrbayes 3.1.2 (Ronquist & Huelsenbeck, 2003). Ten million generations of metropolis-coupled Markov chain Monte Carlo chains with default heating values were used, and the first 50% of sampled trees were discarded as burn-in. Effective sample sizes (ESS) were checked using tracer 1.6 (Rambaut, Suchard, Xie, & Drummond, 2014), and parameters with ESSs lower than 200 were treated as unreliable. We conducted ML analyses in raxMl gUI (Silvestro & Michalak, 2012) with 1,000 bootstrap pseudoreplicates. The partition strategy was as same as BI, and the GTRGAMMA model was implemented for each data partition. paUp* 4.0b10a (Swofford, 2003) was used to perform MP analyses. All characters were treated as unordered and equally weighted. Heuristic searches with tree bisection reconnection were executed for 1,000 random addition replicates. Nodal reliability was assessed using 1,000 bootstrap pseudoreplicates.
Based on the combined mtDNA and nuDNA data, we estimate relationships among different groups using *beast (Heled & Drummond, 2010). RAG-1 was not used as all variations were in heterozygous states. We defined three groups: W, S, and N. Groups W and S were defined based on mtDNA genealogy and distribution. There are three matrilines in nEQTP. Matriline N3 was not included as it was represented by only one individual and we failed to amplify GCG of this sample, although we tried several times. We tested whether samples from nEQTP (matrilines N1 and N2) deviated from Hardy-Weinberg equilibrium (HWE) using arleqUIn 3.5 (Excoffier, Laval, & Schneider, 2005).
If we cannot detect significant deviations, N1 and N2 would be combined as one group (N) in the analyses. Sequences from Quasipaa boulengeri were obtained from another study  and used as the out-group. No shared mtDNA haplotype was found among the three groups. The results of dispersal model test supported that they were completely isolated from each other. Thus, significant violation of the assumption of the method, no gene flow between groups, was not detected in our datasets. A strict molecular clock model was assumed for all loci, and each locus was specified with its own model. We also performed same analyses using nuDNA data only, in case of discordances between mtDNA and nuDNA data. We ran the MCMC chain for 100 million iterations; sampling every 10,000 generation, first 50% trees were discarded as burn-in. Effective sample sizes were checked using tracer 1.6 (Rambaut et al., 2014) to ensure the convergence.
To visualize overall similarity of the haplotypes, we built a medianjoining network (MJN) using netwOrk 4.5.13 (Bandelt, Forster, & Rohl, 1999) based on mitochondrial data. Nuclear data were not involved in this analysis due to a dearth of variation. The MP option was used to remove excessive links and median vectors (Polzin & Daneshmand, 2003).

| Divergence times
Strict molecular clock hypothesis was tested using the likelihood ratio test. We used pHylIp 3.69 (Felsenstein, 2004) to assess the likelihood of unconstrained and clock-enforced matrilineal genealogies. The chisquare test was used to test the significance (Felsenstein, 1981). A priori, we decided to apply this model if the strict molecular clock hypothesis could not be rejected (Ho & Duchêne, 2014).
Times of divergence were estimated in beast 1.8.1 (Drummond, Suchard, Xie, & Rambaut, 2012) based on mitochondrial data. The divergence time between N. pleskei and N. parkeri (8.9 ± 2.7 Ma) (Che et al., 2010) was used as a secondary calibration point. We performed analyses for 50 million generations and sampled trees every 1,000 generations. Burn-in and convergence of the chains were determined with tracer 1.6. The ESS of each parameter was also checked using tracer 1.6 to confirm acceptable convergence.

| Genetic structure and demographic history
Populations of N. pleskei were grouped by using spatial analysis of molecular variation (SAMOVA) (Dupanloup, Schneider, & Excoffier, 2002), which was implemented in spads 1.0 (Dellicour & Mardulyn, 2014). The number of groups (K), which was allowed to range from 2 to 10, was explored with 10,000 iterations and 100 repetitions. The optimum value of K was identified by exploring the behavior of the proportion of total genetic variance due to differences between groups of populations (FCT). We explored multiple K values after FCT reached the plateau (Meirmans, 2015). We also explored the spatial pattern of nuDNA alleles to check whether they were concordant with the population grouping based on mtDNA. To visualize the spatial pattern of nuDNA alleles, we mapped allele frequencies on each locality for three loci separately.
We assessed the demographic history based on mtDNA data of each matrilineal lineage separately, because population subdivision may have influenced detection of demographic shifts. Firstly, we calculated Tajima's D (Tajima, 1989) and Fu's Fs statistics (Fu, 1997) using arleqUIn 3.5 (Excoffier et al., 2005). The significance of each test was assessed by generating null distributions from 10,000 coalescent sim- ulations. Secondly, we tested for signals of demographic expansion by mismatch distributions (Rogers & Harpending, 1992) for each lineage as implemented in arleqUIn 3.5. The sudden demographic expansion was used as the null model, and the null distribution was obtained by generating 10,000 permutations. The significance of deviations from this model was tested by using the sum of squared deviation and raggedness index (Rag). Bioclimatic data for LGM were disaggregated into same resolution as current data. All of the 19 bioclimatic layers were used in analyses because of the controversy about whether correlated variables should be removed or not (Merow, Smith, & Silander, 2013).

| Species distribution model building
We used random null distributions to test whether the SDM of N. pleskei was significantly better than the random distribution (Smith, 2013). Twenty-one random points (the same number as sampling localities) across the study range were used to generate the SDM. We repeated this process 100 times and extracted the values of areas under the curves (AUCs) as the null distribution.

| Testing the barrier effects of rivers
We tested for potential barrier effects of rivers using two approaches.
Firstly, we tested the hypothesis that the rivers did not stop dispersal. For the test, we divided populations into six regions according to the four major rivers in the EQTP ( Figure S1) and calculated F st matrix. Then, we set 17 dispersal models among these groups. Details of each model are summarized in Table 1. Dispersals between populations were estimated by using MIgrate-n 3.2.6 (Beerli, 2006;Beerli & Felsenstein, 2001). Model comparison was done using Bayes factor test (BFT) following Beerli's methods (Beerli & Palczewski, 2010).
The natural log Bayes factors (LBF) were calculated following the formula: LBF = 2 (ln(Prob(D |Model A) -ln(Prob(D|Model i)), while Model A stood for the model with largest log likelihood. LBF larger than 2 suggests Model A was significantly better than Modle i. MIgrate-n was run with one long MCMC chain for 5 million generations while sampling every 1,000 generations and discarding the first 10,000 genealogies as burn-in. Static heating was employed, and temperatures of each Markov chain were set as 1, 1.5, 3, and 10 million. Performance of each run was checked by assessing effective sample size.
We also tested the hypothesis that climate along rivers was suit-  (Hijmans et al., 2012). We checked whether the data vi-

| Matrilineal genealogy and population relationships inference
The best partition strategy of the concatenated mtDNA data was  similar topologies for the matrilineal genealogy (Figures 1 and S2).
Incongruent nodes received weak supports. For example, in MP tree, Matriline S was recovered as a monophyly but not strongly supported.
In BI and ML trees, this monophyly was not recovered. samples from localities 5 and 6 clustered together as did samples from localities 12 and 13. The Jinsha River separated these two sublineages.
Matrilines N1 and W clustered together with strong support.

Significant deviation from Hardy-Weinberg equilibrium (HWE)
was not detected in Tyr (p = .53) and GCG (p = .33). Samples from nEQTP (matrilines N1 and N2) were combined together. The population tree based on the combined mtDNA and nuDNA data ( Figure 2) consisted of mtDNA genealogy (Figure 1 and S2). Population W and N clustered together, and Population S was placed at a basal position.
Population tree based on two nuDNA loci yielded the same topology with tree based on combined data, but with much lower posterior probability ( Figure S3).
The MJN based on the mtDNA ( Figure S4) was congruent with the mtDNA genealogy (Figures 1 and S2). The five major matrilines were depicted. Matriline N1 contained 17 haplotypes and formed a star-like pattern with most common H3 in the center. Haplotypes of Matriline S clustered together.  (Figure 1).

| Genetic structure and demographic history
Population grouping based on mtDNA suggested four groups as FCT

| Species distribution reconstruction
The SDM for N. pleskei differed significantly from the random dis-  (Table 2).  Bioclimatic variables 2, 12, 13, 14, 16, 17, 18, and 19 were deviated from normal distribution. The PCA reduced the bioclimatic variables to three PCs. PC1, PC2, and PC3 explained 58.798%, 23.058%, and 9.176% of the total variance, respectively. The additional components were excluded from subsequent analyses. Scatter plots of these PCs supported that localities of N. pleskei differed from extracted points along the rivers. However, they were mixed with parts of the Yellow and Dadu rivers (Figure 7).

| Different influences of rivers on Nanorana pleskei
Our results reject the null hypotheses that different geographic features in northern EQTP and southern EQTP influenced N. pleskei in a similar way. Our analyses demonstrate clearly that different rivers, or different parts of them, differentially influence N. pleskei. In the sEQTP, rivers tend to be strong barriers to dispersal for N. pleskei, but Genetic structures recovered by matrilineal genealogy are different between nEQTP and sEQTP. Three matrilines, including N1, N2, and N3, are found in nEQTP (Figure 1). But the lineages divergences are not corresponding to the rivers. Shared mtDNA haplotypes are detected in localities that separated by the Yellow River and upstream of Dadu river (Figure 1, Table S1). But in sEQTP, matrilineal genealogy divergence pattern consists of rivers. Matriline S is separated from other matrilines by Dadu and Yalong rivers. Matriline W is separated from its sister lineage, Matriline N1, by Yalong River. Two sublineages are recovered within Matriline W, and they are separated by the Jinsha River ( Figure 1). F st matrix also supports the differences between nEQTP and sEQTP, as F st values between regions separated by rivers in nEQTP are lower than those in sEQTP (Table 3). The conservative nuDNA loci provide little information, but the spatial patterns of alleles are not in conflict with mtDNA lineages. In RAG-1 gene, two private alleles are detected in matrilines S and W. And in Tyr, allele A2 is found in nEQTP and distributes in two sides of the Dadu River. Locus GCG contains more information and supports that Matriline S is divergent from other matrilines. Geographically widespread alleles occur in all nuDNA loci.
Considering the spatial pattern of mtDNA haplotypes and the physiological characters of amphibian (Beebee, 2005), this pattern probably reflects the persistence of ancestral alleles, not widespread gene flow.
Results from SAMOVA also recover different genetic structure between nEQTP and sEQTP. Two groups, N1 and N2, occur in the nEQTP. They present a pattern not associating with rivers. Group N2 contains two localities (1 and 4), which are separated by the Yellow River ( Figure 3). Haplotypes from different matrilines also occur sympatric in localities 4 and 21 (Figure 1). These results indicate recent or ongoing dispersal and gene flow among them. Dispersal modeling supports this scenario, as the models strongly favored by BFT indicate unrestricted mixing among groups in the nEQTP. Two groups, S and W, are in the sEQTP. Dispersal models indicate complete isolation of groups S and W from other groups by rivers. Rivers also promote isolation within groups S and W. Upon setting K = 5, Group S is divided into two groups separated by the Xianshui River-a branch of Yalong River ( Figure S5). These two subgroups do not share mtDNA haplotypes.
Group W consists of two subgroups, which the Jinsha River separates. instead of marsh, are the most common landscape in the sEQTP (You & Yang, 2013

| Demographic history
Our study firmly rejects the null hypothesis that northern and south-   (Qu, Lei, Zhang, & Lu, 2010). Notwithstanding, another factor may also contribute to this apparent stability. Climate changes in the QTP were much less severe during the LGM when compared to other glacial periods (Owen, Caffee, Finkel, & Seong, 2008;Schäfer et al., 2002;Zheng, Xu, & Shen, 2002). The LGM imposed very little effect on demographic history of many species, such as the plateau pika (Ochotona curzoniae) (Ci et al., 2009), Himalayan hemlock (Tsuga dumosa) (Cun & Wang, 2010), and Tibetan antelope (Pantholops hodgsoni) (Du et al., 2010). Compared to these species, N. pleskei appears to have retained stable populations in heterogeneous environments but experienced demographic changes in homogeneous ones.
This result supports the idea that both climatic and geographic features can impact the evolutionary history of a species.

| Implications for conservation
Nanorana pleskei is considered to be near threatened and decreasing in abundance and number of populations (Wang, Annemarie, Muhammad, & Xie, 2004). A major threat to this species appears to be habitat destruction and degradation caused by overgrazing livestock. Our study discovers significantly different population structure between northern and southern populations of N. pleskei in the EQTP.
This result has important implications for the conservation of the species. Populations in the sEQTP appear to be disproportionally important for conservation. These isolated populations have experienced stable populations. Genetic evidence can serve to design evolutionarily significant units and management units (Moritz, 1994). Our results indicate that there are multiple management units in the southern region.
The sEQTP belongs to the Hengduan Mountains, a biodiversity hot spot (Myers et al., 2000). Barriers to dispersal and less severe climatic changes greatly influence the genetic structure and demographic history of N. pleskei. These factors may also contribute to the high biodiversity in this region, as multiple studies indicate species in this region are isolated by river valleys and retain stable demographics. Our results are not only useful for the conservation of this species, but also for all biota in the EQTP.

| CONCLUSIONS
By revealing the genetic structure and population history of N. pleskei, our study documents the impacts of different environments of the northern EQTP and southern EQTP on organisms. Rivers in northern EQTP and southern EQTP influence N. pleskei in different ways. In the nEQTP, rivers are not barriers but rather provide large patches of continuous habitat that facilitate dispersal, which we detect. However, river valleys and their associated environments in the sEQTP host habitat unsuitable for N. pleskei. The demographic histories of N. pleskei in northern and southern regions differ and correspond with different geographic features. In the nEQTP, N. pleskei experienced demographic expansion, which a homogeneous environment promoted. In contrast, complex geographic features in the sEQTP appear to have generated multiple micro-refugia.

CONFLICT OF INTEREST
None declared.

DATA ACCESSIBILITY
All sequences were uploaded to GenBank. Bioclimate data are downloaded from the WorldClim database. Details regarding individual samples are available in Table S1.