Search tips
Search criteria 


Logo of peerjLatest ArticlesFor AuthorsEditorial BoardPeerJPeerJ
PeerJ. 2017; 5: e3770.
Published online 2017 September 11. doi:  10.7717/peerj.3770
PMCID: PMC5598431

Phylogeographic investigation and ecological niche modelling of the endemic frog species Nanorana pleskei revealed multiple refugia in the eastern Tibetan Plateau

Academic Editor: Marta Riutort


The largest plateau Tibetan Plateau supplied an excellent opportunity to investigate the influence of the Pleistocene events on the high-elevation species. To test for the alternative hypotheses of Pleistocene glacial refugia, we used partial sequences of two mitochondrial genes and one nuclear gene to examine the phylogeographic patterns of the endemic frog species Nanorana pleskei across its known range in the eastern Tibetan Plateau, and conducted species distribution modelling (SDM) to explore changes of its distribution range through current and paleo periods. In all data sets, the species was divided into lineage north occupying open plateau platform and lineage south colonizing the mountainous plateau. The divergence of two major clades was estimated at the early Pleistocene. In mtDNA, lineage north contained northeastern and northwestern sublineages, and lineage south had two overlapping-distributed sublineages. Different lineages possessed distinct demographic characteristics, i.e., subdivision in the northeastern sublineage, historical bottleneck effects and recent expansions in the northwestern sublineage and the southeastern sublineage. SDMs depicted that stable suitable habitats had existed in the upper-middle streams of the Yellow River, Dadu River, Jinsha River and Yalong River. These regions were also recognized as the ancestral areas of different lineages. In conclusion, Nanorana pleskei lineages have probably experienced long-term separations. Stable suitable habitats existing in upper-middle streams of major rivers on the eastern Tibetan Plateau and distinct demographic dynamics of different lineages indicated that the lineages possessed independent evolutionary processes in multiple glacial refugia. The findings verified the profound effects of Pleistocene climatic fluctuations on the plateau endemic species.

Keywords: Phylogeography, Tibetan Plateau, Species distribution modelling, Demography, Genetic structure, Multiple refugia


Pleistocene climate fluctuations are believed to cause many temperate species to change distribution range and evolutionary history (Hewitt, 2000, 2004). Allopatric populations in different regions might be glacially isolated and experienced postglacial expansions (Galbreath, Hafner & Zamudio, 2009; Wang et al., 2013). The isolation level contributed to substantial genetic heterogeneity, and high spatial connectivity might be in favor of gene flow and as well, in these processes, demographic dynamics of lineages presented correspondent tendency (Stewart et al., 2010; Peterson & Ammann, 2013). The models have been extensively investigated in Europe and North America (Hewitt, 2004). In East Asia, the issues have been gradually disclosed with increasing phylogeographic researches (Qu et al., 2010; Liu et al., 2013, 2015), but that is far from well understanding them especially in plateau regions.

Glaciations on Asian plateaus were considered to be asynchronous, even within regional mountain-plateau regions (Owen, Finkel & Caffee, 2002; Owen et al., 2008; Thompson, Thompson & Davis, 2006). As the highest and largest plateau on Earth, the Tibetan Plateau extends several million square kilometers and carries area of an average elevation of 4,500 m above sea level (m.a.s.l). Historical uplifts of the plateau promoted geological and climatic transformations in interior basins and large mountains around it (An et al., 2001). During the Pleistocene, different-periods and different-regions glaciations were indicated by increasing evidences (Shi, 2002; Zheng, Xu & Shen, 2002; Owen et al., 2005). Moreover, this region carries habitats that are nearly unique for Oriental and Palearctic organisms that occur in temperate ecosystems and are inhabited by considerable high-elevation endemic species (Zhang, 2009; Myers et al., 2000).

So the Tibetan Plateau supplied a fine biogeographical condition for investigating the influences of Pleistocene climatic oscillations on the evolutionary history of the high-altitude species. Some investigations have focused on species endemic to plateau and surrounding mountains (Cun & Wang, 2010; Hofmann, 2012; Zhou et al., 2013), but conflicting indications were proposed from patterns of different species. The model of refugia-in-periphery has been suggested by studies on several plant and animal species occupying the Tibetan Plateau and peripheral mountains (Frenzel, Bräuning & Adamczyk, 2003; Zhang et al., 2005; Yang et al., 2008). Recently, many works have indicated that multiple refugia exist in the plateau interior (Liu et al., 2015; Hofmann, 2012; Wang et al., 2010). Therefore, more investigations especially plateau endemic species are needed to reveal the historical scenarios.

Amphibians are considered to be indicators of climate changes because their high sensitivity to environmental changes (Bossuyt & Milinkovitch, 2001; Heinicke, Duellman & Hedges, 2007; Fouquet et al., 2012) and their physiological constraints. A small number of anuran species were endemic to high-elevation Tibetan plateau and among them Nanorana pleskei are distributed from 3,300 to 4,500 m.a.s.l, occupying most Hengduan Mountains and eastern Tibetan Plateau (Fei et al., 2009). The distributional range of the species almost covers origin area of many important river systems, i.e., Yellow River, Dadu River, Minjiang River, Yalong River, Jinsha River and Lancang River, and this region present terrain of open plateau platform in the northwest and mountainous plateau in the southeast (Fig. S1). Obviously, the distributional patterns of the species had been certainly influenced by the Pleistocene climate changes, and it is a striking question whether the species had refugia-in-periphery or multiple-refugia-in-plateau during the glacial epoch.

To test the hypotheses, in this study, we investigated the phylogeographic structure and ecological niche modelling of Nanorana pleskei. We sequenced two mitochondrial genes and one nuclear gene of sampling populations across the distribution range of the species. We aimed to: (1) examine the genetic diversity and genetic structure of the species; (2) explore the characteristics of historical demography; (3) explore whether the species retreated to a narrow refugium at the edge of the plateau or occupied multiple refugia on the plateau during the glacial periods. Integrative implements of phylogeographic analyses and species distribution modelling (SDM) allow us to explore architecture of the evolutionary history and demography dynamics of the plateau endemic species.


Sampling and sequencing

A total of 188 individuals were collected from 15 localities mainly spanning across the known range of Nanorana pleskei on the eastern Tibetan Plateau (Table S1; Fig. 1A). Tissues were preserved in 95% Ethanol at −20 °C until DNA extractions were performed. For the phylogenetic analyses, two individuals of the closely-related species Nanorana ventripunctata were collected, and one individual of congeneric species Nanorana parkeri was also included as an outgroup according to previous phylogenetic analyses (Jiang et al., 2005; Che et al., 2010). The Animal Care and Use Committee of Chengdu Institute of Biology, CAS provided full approval for this purely observational research (Number: CIB2014031010). Field experiments were approved by the Management Office of the Zoige Nature Reserve (project number: ZNR201303006).

Figure 1
Sampling localities, phylogenetic tree and TCS haplotype network for Nanorana pleskei.

Total genomic DNA was extracted using a standard phenol–chloroform extraction procedure (Hillis et al., 1996). Two mithochondrial DNA (mtDNA) fragments, the cytochrome oxidase subunit I (COI) and the NADH dehydrogenase subunit 1 (ND1), were amplified for all samples, and fragments of one single-copy nuclear genes, the exon 1 of rhodopsin (Rhod) gene, were amplified for 30 specimens interspersing in major mtDNA lineages (see “Results”). Primers for amplifying COI and ND1 genes were newly designed using the Primer-BLAST tool on NCBI web with default settings. The complete mitochondrial genome of Nanorana pleskei (GenBank Accession no.: HQ324232.1) was downloaded as template for primer-designing. Primer information was shown in Table S2. Amplifications of mtDNA fragments were performed under the following conditions: 95 °C for 3 min; 35 cycles of 95 °C for 35 s, 53 °C (for ND1)/58 °C (for COI) for 35 s, and 72 °C for 70 s; and 72 °C for 10 min. Amplification of the Rhod gene was performed according to the procedures and using primers in the previous studies (Che et al., 2009). PCR products were sequenced in an ABI 3730 sequencer in both directions.

Mithochondrial DNA sequences were assembled in SEQMAN v7.21 (DNAStar Inc., Madison, WI, USA) and aligned using CLUSTAL X v1.83 (Thompson, Gibson & Plewnia, 1997) with default settings and manually checked. Two mtDNA fragments were concatenated and haplotypes were determined using DNASP v5.10.1 (Librado & Rozas, 2009). For the Rhod gene, recombination detection was conducted using seven methods in the program RDP v4.56 (Martin et al., 2015). No signal was found for recombination. Haplotypic state of the Rhod gene was inferred under a Bayesian framework in PHASE (Stephens & Donnelly, 2003) with a probability threshold of 90% and five independent runs. The input files for PHASE were produced by the web program SEQPHASE ( All sequences were submitted to GenBank with accession numbers KX806663, KX806664KX806851, KX806852KX807039, KX807040KX807069, KX807070, KX807071KX807072.

Phylogenetic inferences and dating estimations

Because there were few variable sites in nuclear gene dataset (only three haplotypes were identified in all 30 sequences of the Rhod gene; see “Results”), phylogenetic analyses were conducted only based on mtDNA dataset. MtDNA data were analyzed using maximum likelihood (ML) and Bayesian inference (BI) methods, as implemented in the program PHYML 3.0 (Guindon & Gascuel, 2003) and MRBAYES v3.2.5 (Ronquist et al., 2012), respectively. Before ML and BI analyses, we divided each data set into six partitions through defining the first, second and third codon positions of protein-coding genes, and then used PARTITIONFINDER v1.1.1 (Lanfear et al., 2012) to select best-fitting partitioning schemes and the best-fitting nucleotide substitution model of each partition under the Bayesian information criterion (Schwarz, 1978). The PARTITIONFINDER analysis selected the K80 model for (COI_pos1 + ND1_pos1), F81+I model for (COI_pos2 + ND1_pos2) and TrN model for (COI_pos3 + ND1_pos3). ML and BI analyses were conducted under the selected best-fitting partitioning schemes and the best-fitting nucleotide substitution model of each partition. For ML analyses, the default tree search approach using simultaneous nearest neighbor interchange method and BioNJ tree as starting tree was used to estimate ML tree topologies. Non-parametric bootstrapping with heuristic searches of 1,000 replicates was used to assess confidences of branches in ML trees (Felsenstein, 1985; Felsenstein & Kishino, 1993; Huelsenbeck & Hillis, 1993). For BI analyses, we unlinked parameters for each partition and allowed branch lengths to vary proportionately across partitions. Two independent runs were initiated each with four simultaneous Markov Chain Monte Carlo (MCMC) chains for 30 million generations and sampled every 1,000 generations. Convergence of runs and burn-in period (first 25%) was determined using the program TRACER v1.6 (Rambaut et al., 2014a). A final majority-rule BI tree and the posterior probabilities were achieved from the remaining trees.

Based on mtDNA data, the time to the most recent common ancestor (TMRCA) of Nanorana pleskei was estimated using a Bayesian MCMC approach under a “relaxed molecular clock” model in BEAST v1.8.1 (Drummond et al., 2012). There was no fossil record of Nanorana pleskei that could be used for calibrations. The mutation rate of mitochondrial genes has been found to be broadly constant at 0.57–0.96% change per lineage per million years across many amphibian groups, such as hynobiid salamanders (Weisrock et al., 2001), Bufo (Macey et al., 1998), Ranid frogs (Macey et al., 2001) and Eleutherodactylus toads (Crawford, 2003). The Bufo species were indicated to have a broadly universal substitution rate of 0.65% change per lineage per million years on the mitochondrial ND1–ND2 gene region (Macey et al., 1998), while the Rana boylii species group was also suggested to have a broadly universal substitution rate of 0.65% change per lineage per million years on the mitochondrial ND1, ND2 and CO1 genes (Macey et al., 2001). Accordingly, here we used this mean rate to estimate divergence time among Nanorana pleskei lineages. In BEAST analyses, no partition scheme was selected because we used one mean rate for whole mtDNA data. The GTR + G + I nucleotide substitution model, the most flexible model available in beast, was used to allow the sample space of the parameters to be fully exploited (Huelsenbeck & Rannala, 2004). Genealogy was reconstructed with an uncorrelated lognormal tree prior with a constant population size assumption. A lognormal mean substitution rate of 0.65% change per lineage per million years was used with a lognormal standard deviation of 0.7 producing a 95% credible sampling interval (CI) from 0.5% to 1% change per lineage per million years. The MCMC chains were run for 50 million generations with sampling every 1,000 generations and 10% of the initial samples were discarded as burn-in. The convergence of chains was verified using Tracer by checking that the sampling achieved stationarity and the effective sample size (ESS) for parameters sampled from the MCMC analyses was more than 200. The remaining trees were used to obtain the subsequent maximum clade credibility summary tree with posterior probabilities using TREEANNOTATOR (Rambaut et al., 2014b).

Population genetic structure and demography analyses

Haplotype diversity (Hd) and nucleotide diversity (π) for mtDNA lineages and populations were estimated using DNASP. The significance was tested using 1,000 computer permutations. To visualize geographical patterns of genetic diversity, Hd and π were spatially interpolated using the Kriging method, implemented in the “Geostatistical Analyst” of ARCGIS v10.2 (ESRI, Redlands, CA, USA).

Because the nuclear sequences contain few variable sites and include only a subset of specimens, the relationships between haplotypes of the nuclear gene were only shown using haplotype network. Haplotype networks of mtDNA and the Rhod genes were constructed using maximum parsimony method in the software TCS v1.21 (Clement, Posada & Crandall, 2000).

To access the extent to which the phylo-groups and/or geo-groups explain variation, hierarchical analyses of molecular variance (AMOVA) was conducted in ARLEQUIN v3.1 (Excoffier, Laval & Schneider, 2005). Three kinds of population arrangements were tested: two major lineages, four sublineages and three geographical groups (see “Results”). The significance of the fixation indices Fct, Fsc and Fst derived from AMOVA analyses was determined by 1,000 permutation replicates.

Pairwise Fst between lineages and populations was calculated in ARLEQUIN. To test whether genetic differentiation between populations was derived from the isolation by distance (IBD) model with a significant correlation between geographic and genetic distance (Fst), Mantel tests with 10,000 randomizations were performed in the program IBDWS v3.23 (Jensen, Bohonak & Kelley, 2005).

Genetic signals of departure from neutrality (potential population expansion) were estimated using Tajima’s D (Tajima, 1989) and Fu’ Fs (Fu, 1997) statistics. The values were calculated in Arlequin with 1,000 computer permutations. Mismatch distributions comparing observed and expected distributions of pairwise nucleotide differences were also simulated to test a model of a sudden population expansion. The significance of deviation from this model was evaluated using sum of squares deviations (SSD) and Harpending’s raggedness index (HRI) with significant P values indicating rejection of the recent expansion hypothesis. The analyses were implemented in ARLEQUIN and the significance of deviation of values was tested with 10,000 permutations. The coalescent-based Bayesian skyline plot (BSP), as implemented in BEAST, was used to exhibit demographical fluctuations. The analyses were only conducted for total population and two major lineages because coalescent simulations of sublineages could not reach convergence due to much less variations and samples. In BSP analyses, the starting tree was randomly generated and Bayesian skyline was used for the tree prior. According to the number of samples, two to four grouped coalescent intervals were prior chosen for lineages. The lognormal relaxed clock model and the mutation rate of 0.65%/Ma (95% CI’s: 0.5–1.0%/Ma) as suggested above were specified. MCMC chains were run for 100 million iterations with sampling every 1,000 iterations and discarding the first 10% as burn-in. TRACER was used to check for convergence of chains and ESS.

Biogeographical analyses

The SDM was inferred using maximum entropy machine-learning algorithm in MAXENT v3.3.3k (Phillips, Anderson & Schapire, 2006). Nineteen bioclimatic layers (BIO1-19) describing temperature and precipitation were downloaded from the WorldClim v1.4 database ( (Hijmans et al., 2005). The layers of the last glacial maximum (LGM, ~21,000 years before present) scenarios were obtained by the community climate system model (CCSM; 2.5 arc-minutes resolution; Collins et al., 2006) and the model for interdisciplinary research on climate (MIROC; 2.5 arc-minutes resolution; Hasumi & Emori, 2004). The layers of the current (1950–2000) and last interglacial (LIG; 140,000–120,000 years before present) conditions were retrieved with a resolution of 30 arc-second (Otto-Bliesner et al., 2006), and then were resampled into 2.5 arc-minutes resolution. To eliminate the effect of model over-fitting, Pearson’s correlation coefficient (r; Pearson, 1895) were used to evaluate pairwise correlations between bioclimatic variables. When the absolute value of r > 0.9 indicating that two variables were highly correlated, the one that was more biologically meaningful for the species was chosen for analyses. The results suggested that seven variables (BIO1-4, 12, 14 and 15; see could be chosen for SDMs. Correlation structure among these bioclimatic variables appeared to be stable and consistent across time periods (Mantel test: P < 0.01), suggesting that the model may be transferable to past-time bioclimatic conditions (Jiménez-Valverde et al., 2009).

A total of 55 occurrence records of Nanorana pleskei were selected from the GBIF web site ( and our field records (Table S3). Spatially auto-correlated occurrence points could make models biases in SDM (Veloz, 2009; Hijmans, 2012; Boria et al., 2014). So our occurrence localities were spatially rarefied at 20 km2 based on topographic and climate heterogeneity (Fig. S1), resulting in 46 points for SDM (Table S3). Additionally, a latitudinal background selection bias file accounting for bias associated with latitudinal changes in the studied region and limiting selection of background points to the regions only feasibly colonized by the species was created for SDM (Brown, 2014). In this analysis, the distance to a buffer minimum convex polygon (MCP) (Kremen et al., 2008) was set as two decimal degrees.

To optimize the model performance (Shcheglovitova & Anderson, 2013), different combinations of the five model feature class types (1. linear, 2. linear and quadratic, 3. hinge, 4. linear, quadratic and hinge and 5. linear, quadratic, hinge, product and threshold, as in (Brown, 2014) and a range of regularization parameters (0.5–5; in 0.5 increments)) were tested. The distribution models were calibrated using spatial jackknifing (k = 3) (Brown, 2014), and the best model was chosen with the lowest omission rate and the highest AUC (the area under the receiver operating characteristic curve) (Fielding & Bell, 1997). AUC > 0.8 means that the model performed well (Gassó et al., 2012). The cross-validation strategy with 90% of the presence data as training data and the remaining 10% as test data was used to access the accuracy of models. The model used 10,000 maximum iterations and 10 replications. The 10th percentile training presence logistic threshold commonly used for models comparison (Pearson, 2007) was specified.

The SDMs of all periods were summed and reclassified for predicting stable niche regions (Chan, Brown & Yoder, 2011). To correct over-predictions (Brown, 2014), the final SDMs were clipped by a buffered MCP as above. To visualize the spatial connectivity, dispersal networks among populations were constructed based on the shared mtDNA haplotypes using the weighted least-cost-corridors method in Chan, Brown & Yoder (2011). A friction layer used for least-cost-corridors calculations was created by inverting the niche stable map (Brown, 2014).

The predictions under current and historical bioclimatic conditions were created by MAXENT, and the other operations associated with SDM were implemented in the program SDMTOOLBOX v1.1b (Brown, 2014) applied in ARCGIS.

Ancestral locations were inferred using the BI discrete areas in the software RASP v3.2 (Yu et al., 2015). MCMC trees and the final time tree resulted from the dating estimations above were used as input genealogies. Our 15 sampling locations with GPS coordinates were coded as independent characters. MCMC chains were run for five million iterations with sampling every 1,000 iterations and discarding the first 10% as burn-in. Ancestral sites with their origin ages were interpolated into a continuous map in ArcGIS.


Phylogenetics and dating

The two-gene mtDNA dataset was with length of 1,548 bps, contained 55 variable nucleotide positions, with 44 parsimony informative sites, respectively. ML and BI analyses yielded almost consistent topology (Fig. 1B) except some tips with low supports. ML bootstrap proportions (bp) and Bayesian posterior probabilities (bpp) were more than 95% supporting the monophyly of Nanorana pleskei. In Nanorana pleskei, two major lineages (north and south) were revealed. Lineage north (bp = 71/bpp = 0.96) was comprised of two sublineages, e.g., sublineages northwestern (NW; bp = 94/bpp = 1.00) and northeastern (NE; bp = 50/bpp = 0.94), which were allopatric (Fig. 1A). Lineage south (bp = 88/bpp = 1.00) also contained two sublineages, sublineages S1 (bp = 60/bpp = 0.52) and S2 (bp = 65/bpp = 0.95), which were overlapped in two populations (population 1 and 2; see Fig. 1A).

Mean pairwise uncorrected p-distances within lineages north and south were 0.4% and 0.2%, respectively, much lower than that between them (2.1%). Similarly, the p-distance between sublineages NW and NE was 0.6%, larger than that within them (NW: 0.1%, NE: 0.2%). But the p-distance between sublineages S1 and S2 was 0.2%, little larger than that within them (S1: 0.1%, S2: 0.1%).

The divergence between two major lineages north and south was traced back to about 1.41 million years ago (Mya) with 95% highest posterior density (95% HPD: 2.39–0.85 Mya), and divergence time was 0.55 Mya (95% HPD: 0.78–0.21 Mya) within lineage North and 0.25 Mya (95% HPD: 0.49–0.1 Mya) within lineage south, respectively (Fig. 1B).

Population genetic structure and demographic history

Lineage north and sublineage NE had much higher Hd and π than the sublineage NW, lineage South and sublineages S1 and S2 (Table 1; Fig. 2). As noted, Hd of sublineages S2 was quite low (Table 1).

Table 1
Genetic diversity, neutrality tests and mismatch goodness-of-fit tests in the mtDNA lineages.
Figure 2
Surface of the interpolated genetic diversity of Nanorana pleskei.

TCS analysis of mtDNA data yielded two allopatric and unconnected haplotype groups (95% confidence interval; Fig. 1C), corresponding to the lineages north and south (Fig. 1B). In lineage north, haplotypes from all northeastern populations were mainly clustered into three loops, and haplotypes from all the northwestern populations showed a star-like shape and connected with the northeastern populations in six steps. In lineage south, sublineage S1 had two dominant haplotypes but S2 showed star-like shape. Although the nuclear Rhod gene only had three haplotypes and did not reveal well-supported haplogroups, the northern populations had their own distinct haplotypes as well as the southern populations (Fig. 1D).

Analyses of molecular variance showed significant levels (P < 0.01) across all hierarchical levels when all grouping allocations were tested (Table 2), rejecting the model of no population genetic structure. Highest Fct (95.951%) was found when using grouping arrangement with four sublineages, but considerable Fct was also obtained for grouping arrangements according to two major lineages (Fct = 94.317%) and three geographic groups (Fct = 88.647%). IBDWS analyses resulted in an uncorrelated relationship (r = 0.2593; P = 0.373) between genetic distance (Fst/1 − Fst) and geological distance, rejecting the IBD model.

Table 2
Results of analysis of molecular variance (AMOVA) of mtDNA data.

Tajima’s D value in Neutrality tests was significantly negative (P < 0.05) only in sublineage S2, suggesting past population expansion and significantly negative (P < 0.01) Fu’s FS values were found in sublineages NW and S2 (Table 1). In mismatch distribution, SSD and HRI values suggested population expansion model for the sublineages NE, NW and S2. sublineages NW and S2 exhibited an observed unimodal mismatch frequency distribution (Fig. 3), fitting a recent sudden population expansion model. The BSP for sublineages NW, S1 and S2 was not constructed because of few variations accounting for short coalescent time in each of them. BSPs (Fig. S2) depicted a model of past increase of effective population size in total population of the species (from about 0.05 Mya), lineage north (from about 0.05 Mya) and sublineage NE (from about 0.02 Mya), but this model was not revealed in lineage south suffering short coalescent time with stable population size.

Figure 3
Mismatch distributions of mtDNA lineages.

Species distribution modeling and ancestral locations

The model performances were good and robust because the mean ± standard deviation of AUC was 0.880 ± 0.007 and 0.876 ± 0.048 for the training data and test data, respectively. The SDM under present bioclimatic conditions showed substantial habitat suitability across the current known range of the amphibian species on the eastern Tibetan Plateau (Fig. 4A). The distribution predicted under MIRCO of LGM was similar to the CCSM model (Fig. 4B). Compare to the current distribution, during the LGM, the predicted distribution exhibited obvious contraction to the northeastern region, but in the northwestern region, suitable habitats still existed in the main valley of the upper-middle reaches of Jinsha and Yalong rivers (Fig. 4B). During the LIG, the predictions showed a slight south-direction shift of the distribution range (Fig. 4C). The stable suitable habitats were constructed by summing all SDMs of different periods, and were showed using the highest 30% reclassified values. The results suggested that the stable suitable habitats existed in the regions of eastern half current range of the species (Fig. 4D).

Figure 4
Species distribution models, hypotheses of habitats stability and spatial genetic patterns of Nanorana pleskei.

Because there were no haplotypes between three major geo-groups, spatial connectivity showed by dispersal networks between populations was constructed only within each of them (Fig. 4E). In the southern group, high connectivity was presented near the Yalong River, and in the northeastern group, the highest connectivity existed between the populations near the upper reach of Yellow River, but in the two northwestern populations, the level of connectivity was relatively low.

Ancestral geographical distributions for major mtDNA lineages were inferred to be located in the three geo-groups same as the dispersal-network analyses: region near the middle Yalong River, northeastern region near the upper Yellow River and Dadu River and northwestern region near the upper Yalong and Jinsha valleys (Fig. 4F).


Integrative application of phylogeographic investigation and ecological niche modelling has become a representative approach for explaining intraspecific evolutionary patterns associated with climate oscillations (Chan, Brown & Yoder, 2011; Waltari et al., 2007; Edwards, Keogh & Knowles, 2012). As a case, our study on endemic frog species Nanorana pleskei in the eastern Tibetan Plateau revealed concordance of genetic and ecological modelling signals: first, times of diversification of its lineages occurred through the middle-late Pleistocene, basically fitting scenarios of regionally paleo-climatic fluctuations on the Tibetan Plateau; second, deeply independent origin of major lineages and the spatial interpolations indicated multiple refugia for the species existing across the eastern Tibetan Plateau, also revealed by the niche stability of paleo and current periods across its main distribution range. The integrations and fitting-models allocated substantial insights on evolutionary history of the species and its response to climatic fluctuations.

In Nanorana pleskei, two major mtDNA lineages (north and south) were revealed, showing no shared haplotypes and deep allopatric divergence (ca. 1.41 Mya; Fig. 1B), a similar find as supported by the analyses on the nuclear Rhod gene, i.e., north and south populations had their own haplotypes (Fig. 1D). Recent uplifts of the Tibetan Plateau and Pleistocene glaciations could contribute to the isolations in modern species (Wang et al., 2010; Liu et al., 2015; Zhou et al., 2014). Accordingly, in the range of the lineage North, from ca. 1.2 Mya, the upstream of the Yellow River system had been dramatically upraised (named as the Kunlun-Yellow River tectonic movements; Cui, Wu & Liu, 1997), increasing the landscape differentiations between this region exhibiting open and broad prairie and the range of the lineage south carrying high topographical heterogeneity with tableland and alpine valleys at the southeastern margin of plateau (Fig. 1A; Fig. S1). The movement also promoted expansions of ice-sheets, e.g., the oldest glaciation in the Tibetan Plateau was reported beginning from ca. 1.17 Mya (Shi, Zheng & Su, 2006). Obviously, the early Pleistocene geological and climatic fluctuations would contribute to the isolation and deep divergence of the two lineages. This pattern has also been found in other organisms in this region (Wu et al., 2013; Guo, Liu & Wang, 2012). The divergence of sublineages NW and NE (ca. 0.55 Mya) was likely derived from similar historical factors. Continuous extension of the Kunlun-Yellow River tectonic movement during the middle Pleistocene had promoted the differentiation between hinterland of Tibetan Plateau and edge regions. The range of northwestern lineage (NW) was west of the 600 mm annual precipitation, and adjoined platform of Tibetan Plateau, being more arid and drier, in favor of maintaining its genetic independence from eastern regions. Of course, the maximum glaciation occurred during ca. 0.77–0.5 Mya (Shi, Zheng & Su, 2006) might also provide climate-driving force on separation of it from eastern lineages.

In this view, the species had been likely derived from multiple glacial refugia rather than single refugium. What is more, several lines of evidence support the scenario of multiple glacial refugia at least on the eastern plateau. Firstly, the single-refugium hypothesis generally advocated star-like haplotype network in the whole population and presented decreasing trends on genetic diversity from one single refugium to the newly colonized routes and/or fringe areas (Provan & Bennett, 2008). The pattern was not revealed in Nanorana pleskei. Not only mtDNA data but also nuclear DNA data presented two major distinct haplotype groups (Figs. 1C and and1D),1D), especially in mtDNA, the two lineages independently possessed high genetic diversity (Table 1). As mentioned above, these two lineages dated back to the early Pleistocene much earlier than the LGM were allopatric, implying that they had probably occupied independent glacial refuges. Secondly, the regions with stable suitable environments through history were often inferred as the glacial refugia (Chan, Brown & Yoder, 2011). Our SDMs of current and paleo periods (LIG and LGM) indicated that stable suitable habitats through times for the species existed in almost half of its current distributed range mainly in the basins of major rivers (Fig. 4D). Even in the northern and western range, major valleys, e.g., upper-middle reaches of Yellow River, Dadu River, Yalong River and Jinsha River, still hosted considerable niche stability for the frog (Fig. 4D). Ancestral area distribution also inferred that the MRCA of the populations in these areas were relative old (Fig. 4F). Thus, it is therefore plausible to postulate that these regions had been occupied as the refugia for different lineages. Some other plateau species have likely also taken refuge in the regions (Liu et al., 2013; Zhou et al., 2013). Finally, the multiple-glacial-refugia scenario fits well with the multi-regions and multi-scales glaciations hypotheses during the Pleistocene (Owen et al., 2005; Yang et al., 2008). Moreover, substantial proofs indicated that many districts (e.g., Jinsha River, Yalong River and northeastern region) in the range of Nanorana pleskei were free of ice even in mid-Pleistocene glacial periods on Tibetan plateau (Shi, Zheng & Su, 2006).

Distinct refugia environments would generally create different characteristics of demographic dynamics for the separated lineages (Avise, 2000). For lineage south, the neutrality tests and mismatch distributions all did not show signals of a recent expansion model (Table 1), and BSP showed stable population size through short coalescent time (Fig. S2). But its sublineage S2 significantly fitted the recent expansion model (Table 1), and distributed from central of genetic-diversity map to the edge populations (Fig. 2). This indicated that the lineage south has been probably experienced bottleneck effects, but its sublineage S2 has been probably experienced founder effects and a recent expansion (Zink, 2002). A little is similar for the sublineage NW, the relative older ancestor of it with low genetic diversity and a recent sudden expansion model indicated possible historical bottleneck effects in it (Zink, 2002). However, for the lineages NE possessing high genetic diversity, BSPs showed population increase even through the LGM (Fig. S2) but Fu’s Fs and Tajima’s D did not support the model (Table 1). The interpretation for this discordance probably indicated that the populations underwent expansion but then were recently subdivided, subjected to substantial migration, and/or had experienced historical contractions (Ray, Currat & Excoffier, 2003; Castoe, Spencer & Parkinson, 2007; Bertorelle & Slatkin, 1995).

Of course, there are reasons to be cautious about the findings. In this study, only gene trees based on mtDNA were used to estimate the divergence times and the estimated age of differentiation of Nanorana pleskei lineages predated the time frame used for the SDMs. Compared with species tree methods based on multiple loci, the gene tree method tends to overestimate the divergence times (McCormack et al., 2010). Even though it is true for the findings that the history of major lineages was much probably earlier than the LGM, there still were some other models of demographic patterns, for example, populations had probably been experienced contractions and expansions several times due to glacial cycles during Pleistocene. To now, in view of no fossil information and “incomplete” sampling, we cannot appraise the veracity of estimations on divergence times. So in the future work, “complete” sampling and inclusion of species tree method based on multiple loci even genomic data may be the guarantee to verify the models.

The novel phylogeographic pattern provided important conservation implications for the frog species. The long-term allopatric lineages with limited spatial connectivity (Fig. 4E) should be recognized as independent evolutionarily significant units (ESUs) for future protection (Crandall et al., 2000). The drainage basins in the range of the species hosted important water source for the Yellow River and Yangtze River systems and were indicated as the refugia for it. These valleys have been prominently disturbed by human activity especially in the new century. Thus, more effective measures and managements are desired for protecting habitats in the refuges.


This study integrating phylogeographic investigation and ecological niche modelling suggested that Pleistocene climatic fluctuations combined with the uplifts of Tibetan Plateau to profoundly impacted the evolutionary and distributional patterns of the plateau endemic amphibian species. Nanorana pleskei lineages have taken refuge in upper-middle streams of the Yellow River, Dadu River, Jinsha River and Yalong River because the regions have probably loaded stable, long-term suitable habitats. Multiple refugia environments have probably promoted the lineages possessing distinct demographic characteristics. These points provide useful conservation implications for the species.

Supplemental Information


Supplemental Information 1

Genetic diversity and results of Neutrality tests and Mismatch distribution analyses for each MtDNA lineages.

Supplemental Information 2

Primers designed in this study for amplifying and sequencing COI and ND1 genes.

Supplemental Information 3

Occurrence data for species distribution modelling.

Supplemental Information 4

Climatic and topographic heterogeneity in the eastern Tibetan Plateau.

(A) Climatic heterogeneity. (B) Topographic heterogeneity. Warmer color means higher heterogeneity. (C) A landscape photo of the northern region, showing open landscape in prairie. (D) A landscape photo of the southern region, showing narrow landscape in mountains and valleys landscape.


Supplemental Information 5

Bayesian Skyline plots for some lineages of N. pleskei.

The bold black line indicates the median value of effective population size; the thin black lines denote the 95% highest posterior probability interval. The y-axis correspond to population size Ne* Tau (effective population time x generation length time in millions of years).


Supplemental Information 6

Raw data sequences.

Funding Statement

This work was supported by the National Key Research and Development Program of China (2017YFC0505202) and the National Natural Science Foundation of China (Nos. 31201702 and 31471964). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

The following grant information was disclosed by the authors:

National Key Research and Development Program of China 2017YFC0505202.
National Natural Science Foundation of China Nos. 31201702 and 31471964.

Additional Information and Declarations

Competing Interests

The authors declare that they have no competing interests.

Author Contributions

Bin Wang conceived and designed the experiments, performed the experiments, analyzed the data, contributed reagents/materials/analysis tools, wrote the paper and prepared figures and/or tables.

Feng Xie conceived and designed the experiments, wrote the paper and reviewed drafts of the paper.

Jiannan Li performed the experiments, analyzed the data and contributed reagents/materials/analysis tools.

Gang Wang performed the experiments, analyzed the data and contributed reagents/materials/analysis tools.

Cheng Li analyzed the data, contributed reagents/materials/analysis tools.

Jianping Jiang conceived and designed the experiments, wrote the paper, prepared figures and/or tables and reviewed drafts of the paper.

Animal Ethics

The following information was supplied relating to ethical approvals (i.e., approving body and any reference numbers):

The Animal Care and Use Committee of Chengdu Institute of Biology, CAS provided full approval for this purely observational research (CIB2014031010).

Field Study Permissions

The following information was supplied relating to field study approvals (i.e., approving body and any reference numbers):

Field experiments were approved by the Management Office of the Zoige Nature Reserve (project number: ZNR201303006).

DNA Deposition

The following information was supplied regarding the deposition of DNA sequences:

The sequences described here are accessible via GenBank accession numbers: KX806663, KX806664KX806851, KX806852KX807039, KX807040KX807069, KX807070, KX807071KX807072. We have also uploaded the sequences in a Supplemental File.

Data Availability

The following information was supplied regarding data availability:

The raw data has been supplied as Supplemental Dataset Files.


An et al. (2001) An Z, Kutzbach JE, Prell WL, Porter SC. Evolution of Asian monsoons and phased uplift of the Himalaya–Tibetan plateau since Late Miocene times. Nature. 2001;411:62–66. doi: 10.1038/35075035. [PubMed] [Cross Ref]
Avise (2000) Avise JC. Phylogeography: The History and Formation of Species. Cambridge: Harvard University Press; 2000.
Bertorelle & Slatkin (1995) Bertorelle G, Slatkin M. The number of segregating sites in expanding human populations, with implications for estimates of demographic parameters. Molecular Biology and Evolution. 1995;12(5):887–892. doi: 10.1093/oxfordjournals.molbev.a040265. [PubMed] [Cross Ref]
Boria et al. (2014) Boria RA, Olson LE, Goodman SM, Anderson RA. Spatial filtering to reduce sampling bias can improve the performance of ecological niche models. Ecological Modelling. 2014;275:73–77. doi: 10.1016/j.ecolmodel.2013.12.012. [Cross Ref]
Bossuyt & Milinkovitch (2001) Bossuyt F, Milinkovitch MC. Amphibians as indicators of early Tertiary “out-of-India” dispersal of vertebrates. Science. 2001;292(5514):93–95. doi: 10.1126/science.1058875. [PubMed] [Cross Ref]
Brown (2014) Brown JL. SDMtoolbox: a python-based GIS toolkit for landscape genetic, biogeographic and species distribution model analyses. Methods in Ecology Evolution. 2014;5(7):694–700. doi: 10.1111/2041-210x.12200. [Cross Ref]
Castoe, Spencer & Parkinson (2007) Castoe TA, Spencer CL, Parkinson CL. Phylogeographic structure and historical demography of the western diamondback rattlesnake (Crotalus atrox): a perspective on North American desert biogeography. Molecular Phylogenetics and Evolution. 2007;42(1):193–212. doi: 10.1016/j.ympev.2006.07.002. [PubMed] [Cross Ref]
Chan, Brown & Yoder (2011) Chan LM, Brown JL, Yoder AD. Integrating statistical genetic and geospatial methods brings new power to phylogeography. Molecular Phylogenetics and Evolution. 2011;59(2):523–537. doi: 10.1016/j.ympev.2011.01.020. [PubMed] [Cross Ref]
Che et al. (2009) Che J, Hu JS, Zhou WW, Murphy RW, Papenfuss TJ, Chen MY, Rao DQ, Li PP, Zhang YP. Phylogeny of the Asian spiny frog tribe Paini (Family Dicroglossidae) sensu Dubois. Molecular Phylogenetics and Evolution. 2009;50(1):59–73. doi: 10.1016/j.ympev.2008.10.007. [PubMed] [Cross Ref]
Che et al. (2010) Che J, Zhou WW, Hu JS, Yan F, Papenfuss TJ, Wake DB, Zhang YP. Spiny frogs (Paini) illuminate the history of the Himalayan region and Southeast Asia. Proceedings of the National Academy of Sciences of the United States of America. 2010;107(31):13765–13770. doi: 10.1073/pnas.1008415107. [PubMed] [Cross Ref]
Clement, Posada & Crandall (2000) Clement M, Posada D, Crandall KA. TCS: a computer program to estimate gene genealogies. Molecular Ecology. 2000;9(10):1657–1659. doi: 10.1046/j.1365-294x.2000.01020.x. [PubMed] [Cross Ref]
Collins et al. (2006) Collins WD, Rasch PJ, Boville BA, Hack JJ, McCaa JR, Williamson DL, Briegleb BP, Bitz CM, Lin SJ, Zhang MH. The formulation and atmospheric simulation of the Community Atmosphere Model Version 3 (CAM3) Journal of Climate. 2006;19(11):2144–2161. doi: 10.1175/jcli3760.1. [Cross Ref]
Crandall et al. (2000) Crandall KA, Bininda-Emonds ORP, Mace GM, Wayne RK. Considering evolutionary processes in conservation biology. Trends in Ecology & Evolution. 2000;15(7):290–295. doi: 10.1016/s0169-5347(00)01876-0. [PubMed] [Cross Ref]
Crawford (2003) Crawford AJ. Relative rates of nucleotide substitution in frogs. Journal of Molecular Evolution. 2003;57(6):636–641. doi: 10.1007/s00239-003-2513-7. [PubMed] [Cross Ref]
Cui, Wu & Liu (1997) Cui ZJ, Wu YQ, Liu GN. Discovery and character of the Kunlun-Yellow River Movement. Chinese Science Bulletin. 1997;42(18):1986–1989.
Cun & Wang (2010) Cun YZ, Wang XQ. Plant recolonization in the Himalaya from the southeastern Qinghai–Tibetan Plateau: geographical isolation contributed to high population differentiation. Molecular Phylogenetics and Evolution. 2010;56:972–982. doi: 10.1016/j.ympev.2010.05.007. [PubMed] [Cross Ref]
Drummond et al. (2012) Drummond AJ, Suchard MA, Xie D, Rambaut A. Bayesian phylogenetics with BEAUti and the BEAST 1.7. Molecular Biology and Evolution. 2012;29(8):1969–1973. doi: 10.1093/molbev/mss075. [PMC free article] [PubMed] [Cross Ref]
Edwards, Keogh & Knowles (2012) Edwards DL, Keogh JS, Knowles LL. Effects of vicariant barriers, habitat stability, population isolation and environmental features on species divergence in the southwestern Australian coastal reptile community. Molecular Ecology. 2012;21(15):3809–3822. doi: 10.1111/j.1365-294x.2012.05637.x. [PubMed] [Cross Ref]
Excoffier, Laval & Schneider (2005) Excoffier L, Laval G, Schneider S. Arlequin version 3.0: an integrated software package for population genetics data analysis. Evolutionary Bioinformatics. 2005;1:47–50. [PMC free article] [PubMed]
Fei et al. (2009) Fei L, Hu SQ, Ye CY, Huang YZ. Fauna Sinica, Amphibia: Volume 3. Anura: Ranidae. Beijing: Science Press; 2009.
Felsenstein (1985) Felsenstein J. Confidence limits on phylogenies: an approach using the bootstrap. Evolution. 1985;39(4):783–791. doi: 10.2307/2408678. [PubMed] [Cross Ref]
Felsenstein & Kishino (1993) Felsenstein J, Kishino H. Is there something wrong with the bootstrap on phylogeny? A reply to Hillis and Bull. Systematic Biology. 1993;42(2):193–200. doi: 10.2307/2992541. [Cross Ref]
Fielding & Bell (1997) Fielding AH, Bell JF. A review of methods for the assessment of prediction errors in conservation presence/absence models. Environmental Conservation. 1997;24(1):38–49. doi: 10.1017/s0376892997000088. [Cross Ref]
Fouquet et al. (2012) Fouquet A, Noonan BP, Rodrigues MT, Pech N, Gilles A, Gemmell NJ. Multiple quaternary refugia in the eastern Guiana Shield revealed by comparative phylogeography of 12 frog species. Systematic Biology. 2012;61(3):461–489. doi: 10.1093/sysbio/syr130. [PubMed] [Cross Ref]
Frenzel, Bräuning & Adamczyk (2003) Frenzel B, Bräuning A, Adamczyk S. On the problem of possible last-glacial forest-refuge-areas within the deep valleys of eastern Tibet. Erdkunde. 2003;57(3):182–198. doi: 10.3112/erdkunde.2003.03.02. [Cross Ref]
Fu (1997) Fu YX. Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics. 1997;147:915–925. [PubMed]
Galbreath, Hafner & Zamudio (2009) Galbreath KE, Hafner DJ, Zamudio KR. When cold is better: climate-driven elevation shifts yield complex patterns of diversification and demography in an alpine specialist (American pika, Ochotona princeps) Evolution. 2009;63(11):2848–2863. doi: 10.1111/j.1558-5646.2009.00803.x. [PubMed] [Cross Ref]
Gassó et al. (2012) Gassó N, Thuiller W, Pino J, Vilà M. Potential distribution range of invasive plant species in Spain. NeoBiota. 2012;12:25–40. doi: 10.3897/neobiota.12.2341. [Cross Ref]
Guindon & Gascuel (2003) Guindon S, Gascuel O. A simple, fast, and accurate algorithm to estimate large phylogenies by maximum likelihood. Systematic Biology. 2003;52(5):696–704. doi: 10.1080/10635150390235520. [PubMed] [Cross Ref]
Guo, Liu & Wang (2012) Guo XG, Liu L, Wang YZ. Phylogeography of the Phrynocephalus vlangalii species complex in the upper reaches of the Yellow River inferred from mtDNA ND4-tRNALEU segments. Asian Herpetological Research. 2012;3(1):52–68. doi: 10.3724/sp.j.1245.2012.00052. [Cross Ref]
Hasumi & Emori (2004) Hasumi H, Emori S. K-1 Coupled GCM (MIROC) Description. Tokyo: Center for Climate System Research, University of Tokyo; 2004.
Heinicke, Duellman & Hedges (2007) Heinicke MP, Duellman WE, Hedges SB. Major Caribbean and Central American frog faunas originated by ancient oceanic dispersal. Proceedings of the National Academy of the Sciences of the United States of America. 2007;104(24):10092–10097. doi: 10.1073/pnas.0611051104. [PubMed] [Cross Ref]
Hewitt (2000) Hewitt GM. The genetic legacy of the quaternary ice ages. Nature. 2000;405(6789):907–913. [PubMed]
Hewitt (2004) Hewitt GM. Genetic consequences of climatic oscillations in the quaternary. Philosophical Transactions of the Royal Society B: Biological Sciences. 2004;359(1442):183–195. doi: 10.1098/rstb.2003.1388. [PMC free article] [PubMed] [Cross Ref]
Hijmans (2012) Hijmans RJ. Cross-validation of species distribution models: removing spatial sorting bias and calibration with a null model. Ecology. 2012;93(3):679–688. doi: 10.1890/11-0826.1. [PubMed] [Cross Ref]
Hijmans et al. (2005) Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A. Very high resolution interpolated climate surfaces for global land areas. International Journal of Climatology. 2005;25(15):1965–1978. doi: 10.1002/joc.1276. [Cross Ref]
Hillis et al. (1996) Hillis DM, Mable BK, Larson A, Davis SK, Zimmer EA. Nucleic acids IV: sequencing and cloning. In: Hillis D, Moritz C, Mable B, editors. Molecular Systematics. Massachusetts: Sinauer Associates; 1996. pp. 321–378.
Hofmann (2012) Hofmann S. Population genetic structure and geographic differentiation in the hot spring snake Thermophis baileyi (Serpentes, Colubridae): indications for glacial refuges in southern-central Tibet. Molecular Phylogenetics and Evolution. 2012;63(2):396–406. doi: 10.1016/j.ympev.2012.01.014. [PubMed] [Cross Ref]
Huelsenbeck & Hillis (1993) Huelsenbeck JP, Hillis DM. Success of phylogenetic methods in the four-taxon case. Systematic Biology. 1993;42(3):247–264. doi: 10.2307/2992463. [Cross Ref]
Huelsenbeck & Rannala (2004) Huelsenbeck JP, Rannala B. Frequentist properties of Bayesian posterior probabilities of phylogenetic trees under simple and complex substitution models. Systematic Biology. 2004;53(6):904–913. doi: 10.1080/10635150490522629. [PubMed] [Cross Ref]
Jensen, Bohonak & Kelley (2005) Jensen JL, Bohonak AJ, Kelley ST. Isolation by distance, web service v.3.14. BMC Genetics. 2005;6:13. doi: 10.1186/1471-2156-6-13. [PMC free article] [PubMed] [Cross Ref]
Jiang et al. (2005) Jiang JP, Dubois A, Ohler A, Tillier A, Chen XH, Xie F, Stöck M. Phylogenetic relationships of the tribe Paini (Amphibia, Anura, Ranidae) based on partial sequences of mitochondrial 12s and 16s rRNA genes. Zoological Science. 2005;22(3):353–362. doi: 10.2108/zsj.22.353. [PubMed] [Cross Ref]
Jiménez-Valverde et al. (2009) Jiménez-Valverde A, Nakazawa Y, Lira-Noriega A, Peterson AT. Environmental correlation structure and ecological niche model projections. Biodiversity Informatics. 2009;6(1):28–35. doi: 10.17161/bi.v6i1.1634. [Cross Ref]
Kremen et al. (2008) Kremen C, Cameron A, Moilanen A, Phillips SJ, Thomas CD, Beentje H, Dransfield J, Fisher BL, Glaw F, Good TC, Harper GJ, Hijmans RJ, Lees DC, Louis E, Jr, Nussbaum RA, Raxworthy CJ, Razafimpahanana A, Schatz GE, Vences M, Vieites DR, Wright PC, Zjhra ML. Aligning conservation priorities across taxa in Madagascar with high-resolution planning tools. Science. 2008;320(5873):222–226. doi: 10.1126/science.1155193. [PubMed] [Cross Ref]
Lanfear et al. (2012) Lanfear R, Calcott B, Ho SYW, Guindon S. PartitionFinder: combined selection of partitioning schemes and substitution models for phylogenetic analyses. Molecular Biology and Evolution. 2012;29(6):1695–1701. doi: 10.1093/molbev/mss020. [PubMed] [Cross Ref]
Librado & Rozas (2009) Librado P, Rozas J. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics. 2009;25:1451–1452. doi: 10.1093/bioinformatics/btp187. [PubMed] [Cross Ref]
Liu et al. (2013) Liu J, Möller M, Provan J, Gao LM, Poudel RC, Li DZ. Geological and ecological factors drive cryptic speciation of yews in a biodiversity hotspot. New Phytologist. 2013;199(4):1093–1108. doi: 10.1111/nph.12336. [PubMed] [Cross Ref]
Liu et al. (2015) Liu J, Wang CM, Fu DL, Hu XJ, Xie XM, Liu PF, Zhang Q, Li MH. Phylogeography of Nanorana parkeri (Anura: Ranidae) and multiple refugia on the Tibetan Plateau revealed by mitochondrial and nuclear DNA. Scientific Reports. 2015;5(1):9857. doi: 10.1038/srep09857. [PMC free article] [PubMed] [Cross Ref]
Macey et al. (1998) Macey JR, Schulte JA, II, Larson A, Fang Z, Wang Y, Tuniyev BS, Papenfuss TJ. Phylogenetic relationships of toads in the Bufo bufo species group from the eastern escarpment of the Tibetan Plateau: a case of vicariance and dispersal. Molecular Phylogenetics and Evolution. 1998;9(1):80–87. doi: 10.1006/mpev.1997.0440. [PubMed] [Cross Ref]
Macey et al. (2001) Macey JR, Strasburg J, Brisson J, Vredenburg V, Jennings M, Larson A. Molecular phylogenetics of western North American frogs of the Rana boylii species group. Molecular Phylogenetics and Evolution. 2001;19(1):131–143. doi: 10.1006/mpev.2000.0908. [PubMed] [Cross Ref]
Martin et al. (2015) Martin DP, Murrell B, Golden M, Khoosal A, Muhire B. RDP4: detection and analysis of recombination patterns in virus genomes. Virus Evolution. 2015;1(1):1–5. doi: 10.1093/ve/vev003. [PMC free article] [PubMed] [Cross Ref]
McCormack et al. (2010) McCormack JE, Heled J, Delaney KS, Peterson AT, Knowles LL. Calibrating divergence times on species trees versus gene trees: implications for speciation history of Aphelocoma jays. Evolution. 2010;65(1):184–202. doi: 10.1111/j.1558-5646.2010.01097.x. [PubMed] [Cross Ref]
Myers et al. (2000) Myers N, Mittermeier RA, Mittermeier CG, Da Fonseca GA, Kent J. Biodiversity hotspots for conservation priorities. Nature. 2000;403:853–858. doi: 10.1038/35002501. [PubMed] [Cross Ref]
Otto-Bliesner et al. (2006) Otto-Bliesner BL, Marshall SJ, Overpeck JT, Miller GH, Hu A. Simulating arctic climate warmth and icefield retreat in the last interglaciation. Science. 2006;311(5768):1751–1753. doi: 10.1126/science.1120808. [PubMed] [Cross Ref]
Owen et al. (2008) Owen LA, Caffee MW, Finkel RC, Seong YB. Quaternary glaciation of the Himalayan-Tibetan orogen. Journal of Quaternary Science. 2008;23(6–7):513–531. doi: 10.1002/jqs.1203. [Cross Ref]
Owen et al. (2005) Owen LA, Derbyshire E, Finkel RC, Barnard PL, Ma H, Asahi K, Caffee MW, Derbyshireet E. Climatic and topographic controls on the style and timing of late quaternary glaciation throughout Tibet and the Himalaya defined by 10Be cosmogenic radionuclide surface exposure dating. Quaternary Science Review. 2005;24(12–13):1391–1411. doi: 10.1016/j.quascirev.2004.10.014. [Cross Ref]
Owen, Finkel & Caffee (2002) Owen LA, Finkel RC, Caffee MW. A note on the extent of glaciation throughout the Himalaya during the global last glacial maximum. Quaternary Science Review. 2002;21(1–3):147–157. doi: 10.1016/s0277-3791(01)00104-4. [Cross Ref]
Pearson (2007) Pearson RG. American Museum of Natural History; 2007. Species’ distribution modeling for conservation educators and practitioners. Synthesis.
Pearson (1895) Pearson K. Notes on regression and inheritance in the case of two parents. Proceedings of the Royal Society of London (1854–1905) 1895;58(1):240–242. doi: 10.1098/rspl.1895.0041. [Cross Ref]
Peterson & Ammann (2013) Peterson AT, Ammann CM. Global patterns of connectivity and isolation of populations of forest bird species in the late Pleistocene. Global Ecology and Biogeography. 2013;22(5):596–606. doi: 10.1111/geb.12010. [Cross Ref]
Phillips, Anderson & Schapire (2006) Phillips SJ, Anderson RP, Schapire RE. Maximum entropy modeling of species geographic distributions. Ecological Modelling. 2006;190(3–4):231–259. doi: 10.1016/j.ecolmodel.2005.03.026. [Cross Ref]
Provan & Bennett (2008) Provan J, Bennett KD. Phylogeographic insights into cryptic glacial refugia. Trends in Ecology and Evolution. 2008;23(10):564–571. doi: 10.1016/j.tree.2008.06.010. [PubMed] [Cross Ref]
Qu et al. (2010) Qu Y, Lei F, Zhang R, Lu X. Comparative phylogeography of five avian species: implications for Pleistocene evolutionary history in the Qinghai–Tibetan plateau. Molecular Ecology. 2010;19(2):338–351. doi: 10.1111/j.1365-294x.2009.04445.x. [PubMed] [Cross Ref]
Rambaut et al. (2014a) Rambaut A, Suchard MA, Xie D, Drummond AJ. Tracer v1.6 2014a
Rambaut et al. (2014b) Rambaut A, Suchard MA, Xie D, Drummond AJ. Treeannotator 2014b
Ray, Currat & Excoffier (2003) Ray N, Currat M, Excoffier L. Intra-deme molecular diversity in spatially expanding populations. Molecular Biology and Evolution. 2003;20(1):76–86. doi: 10.1093/molbev/msg009. [PubMed] [Cross Ref]
Ronquist et al. (2012) Ronquist F, Teslenko M, van der Mark P, Ayres D, Darling A, Höhna S, Larget B, Liu L, Sichard MA, Huelsenbeck JP. MrBayes 3.2: efficient Bayesian phylogenetic inference and model choice across a large model space. Systematic Biology. 2012;61(3):539–542. doi: 10.1093/sysbio/sys029. [PMC free article] [PubMed] [Cross Ref]
Schwarz (1978) Schwarz GE. Estimating the dimension of a model. Annals of Statistics. 1978;6(2):461–464. doi: 10.1214/aos/1176344136. [Cross Ref]
Shcheglovitova & Anderson (2013) Shcheglovitova M, Anderson RP. Estimating optimal complexity for ecological niche models: a jackknife approach for species with small sample sizes. Ecological Modelling. 2013;269:9–17. doi: 10.1016/j.ecolmodel.2013.08.011. [Cross Ref]
Shi (2002) Shi Y. Characteristics of late Quaternary monsoonal glaciation on the Tibetan Plateau and in East Asia. Quaternary International. 2002;97–98:79–91. doi: 10.1016/s1040-6182(02)00053-8. [Cross Ref]
Shi, Zheng & Su (2006) Shi Y, Zheng B, Su Z. The Quaternary Glaciations and Environmental Variations in China. Shijiazhuang: Hebei Science and Technology Press; 2006.
Stephens & Donnelly (2003) Stephens M, Donnelly P. A comparison of Bayesian methods for haplotype reconstruction from population genotype data. American Journal of Human Genetics. 2003;73(5):1162–1169. doi: 10.1086/379378. [PubMed] [Cross Ref]
Stewart et al. (2010) Stewart JR, Lister AM, Barnes I, Dalén L. Refugia revisited: individualistic responses of species in space and time. Proceedings of the Royal Society B: Biological Sciences. 2010;277(1682):661–671. doi: 10.1098/rspb.2009.1272. [PMC free article] [PubMed] [Cross Ref]
Tajima (1989) Tajima F. Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics. 1989;123(3):585–595. [PubMed]
Thompson, Gibson & Plewnia (1997) Thompson J, Gibson TJ, Plewnia F. The CLUSTAL X windows interface flexible strategies for multiple sequence alignment aided by quality analysis tools. Nucleic Acids Research. 1997;25(25):4876–4882. [PMC free article] [PubMed]
Thompson, Thompson & Davis (2006) Thompson LG, Thompson EM, Davis ME. Ice core evidence for asynchronous glaciation on the Tibetan Plateau. Quaternary International. 2006;154–155:3–10. doi: 10.1016/j.quaint.2006.02.001. [Cross Ref]
Veloz (2009) Veloz SD. Spatially autocorrelated sampling falsely inflates measures of accuracy for presenceonly niche models. Journal of Biogeography. 2009;36(12):2290–2299. doi: 10.1111/j.1365-2699.2009.02174.x. [Cross Ref]
Waltari et al. (2007) Waltari E, Hijmans RJ, Peterson AT, Ny ari AS, Perkins SL, Guralnick RP. Locating Pleistocene refugia: comparing phylogeographic and ecological niche model predictions. PLOS ONE. 2007;2(7):e563 doi: 10.1371/journal.pone.0000563. [PMC free article] [PubMed] [Cross Ref]
Wang et al. (2010) Wang H, Qing L, Sun K, Lu F, Wang YG, Song ZP, Wu QH, Chen JK, Zhang WJ. Phylogeographic structure of Hippophae tibetana (Elaeagnaceae) highlights the highest microrefugia and the rapid uplift of the Qinghai–Tibetan Plateau. Molecular Ecology. 2010;19(14):2964–2979. doi: 10.1111/j.1365-294x.2010.04729.x. [PubMed] [Cross Ref]
Wang et al. (2013) Wang WJ, McKay BD, Dai CY, Zhao N, Zhang RY, Qu YH, Song G, Li S, Liang W, Yang XJ, Pasquet E, Lei FM. Glacial expansion and diversification of an East Asian montane bird, the green-backed tit (Parus monticolus) Journal of Biogeography. 2013;40(6):1156–1169. doi: 10.1111/jbi.12055. [Cross Ref]
Weisrock et al. (2001) Weisrock DW, Macey JR, Ugurtas IH, Larson A, Papenfuss TJ. Molecular phylogenetics and historical biogeography among salamandrids of the “true” salamander clade: rapid branching of numerous highly divergent lineages in Mertensiella luschani associated with the rise of Anatolia. Molecular Phylogenetics and Evolution. 2001;18(3):434–448. doi: 10.1006/mpev.2000.0905. [PubMed] [Cross Ref]
Wu et al. (2013) Wu X, Luo J, Huang S, Chen Z, Xiao H, Zhang Y. Molecular phylogeography and evolutionary history of Poropuntius huangchuchieni (Cyprinidae) in southwest China. PLOS ONE. 2013;8(11):e79975 doi: 10.1371/journal.pone.0079975. [PMC free article] [PubMed] [Cross Ref]
Yang et al. (2008) Yang FS, Li YF, Ding X, Wang XQ. Extensive population expansion of Pedicularis longiflora (Orobanchaceae) on the Qinghai–Tibetan Plateau and its correlation with the quaternary climate change. Molecular Ecology. 2008;17(23):5135–5145. doi: 10.1111/j.1365-294x.2008.03976.x. [PubMed] [Cross Ref]
Yu et al. (2015) Yu Y, Harris AJ, Blair C, He XJ. RASP (reconstruct ancestral state in phylogenies): a tool for historical biogeography. Molecular Phylogenetics and Evolution. 2015;87:46–49. doi: 10.1016/j.ympev.2015.03.008. [PubMed] [Cross Ref]
Zhang (2009) Zhang RZ. Zoogeography of China. Beijing: Science Press; 2009.
Zhang et al. (2005) Zhang Q, Chiang TY, George M, Liu JQ, Abbott RJ. Phylogeography of the Qinghai-Tibetan Plateau endemic Juniperus przewalskii (Cupressaceae) inferred from chloroplast DNA sequence variation. Molecular Ecology. 2005;14(11):3513–3524. doi: 10.1111/j.1365-294x.2005.02677.x. [PubMed] [Cross Ref]
Zheng, Xu & Shen (2002) Zheng B, Xu Q, Shen Y. The relationship between climate change and quaternary glacial cycles on the Qinghai–Tibetan Plateau: review and speculation. Quaternary International. 2002;97–98:93–101. doi: 10.1016/s1040-6182(02)00054-x. [Cross Ref]
Zhou et al. (2013) Zhou WW, Yan F, Fu J, Wu SF, Murphy RW, Che J, Zhang YP. River islands, refugia and genetic structuring in the endemic brown frog Rana kukunoris (Anura, Ranidae) of the Qinghai–Tibetan Plateau. Molecular Ecology. 2013;22(1):130–142. doi: 10.1111/mec.12087. [PubMed] [Cross Ref]
Zhou et al. (2014) Zhou WW, Zhang BL, Chen HM, Jin JQ, Yang JX, Wang YY, Jiang K, Murphy RW, Zhang YP, Che J. DNA barcodes and species distribution models evaluate threats of global climate changes to genetic diversity: a case study from Nanorana parkeri (Anura: Dicroglossidae) PLOS ONE. 2014;9(8):e103899 doi: 10.1371/journal.pone.0103899. [PMC free article] [PubMed] [Cross Ref]
Zink (2002) Zink RM. Methods in comparative phylogeography, and their application to studying evolution in the North American Aridlands. Integrative and Comparative Biology. 2002;42(5):953–959. doi: 10.1093/icb/42.5.953. [PubMed] [Cross Ref]

Articles from PeerJ are provided here courtesy of PeerJ, Inc