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

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).

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 (http://seqphase.mpg.de/seqphase/). All sequences were submitted to GenBank with accession numbers KX806663, KX806664–KX806851, KX806852–KX807039, KX807040–KX807069, KX807070, KX807071–KX807072.

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 (http://www.worldclim.org) (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 www.worldclim.org/bioclim) 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 (www.gbif.es) 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 km² 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).

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.

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.

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).

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).