Figures
Abstract
Background
The influence of Pleistocene climatic fluctuations on intraspecific diversification in the Qinling–Daba Mountains of East Asia remains poorly investigated. We tested hypotheses concerning refugia during the last glacial maximum (LGM) in this region by examining the phylogeography of the swelled vent frog (Feirana quadranus; Dicroglossidae, Anura, Amphibia).
Methodology/Principal Findings
We obtained complete mitochondrial ND2 gene sequences of 224 individuals from 34 populations of Feirana quadranus for phylogeographic analyses. Additionally, we obtained nuclear tyrosinase gene sequences of 68 F. quadranus, one F. kangxianensis and three F. taihangnica samples to test for mitochondrial introgression among them. Phylogenetic analyses based on all genes revealed no introgression among them. Phylogenetic analyses based on ND2 datasets revealed that F. quadranus was comprised of six lineages which were separated by deep valleys; the sole exception is that the Main Qinling and Micang–Western Qinling lineages overlap in distribution. Analyses of population structure indicated restricted gene flow among lineages. Coalescent simulations and divergence dating indicated that the basal diversification within F. quadranus may be associated with the dramatic uplifts of the Tibetan Plateau during the Pliocene. Coalescent simulations indicated that Wuling, Daba, and Western Qinling–Micang–Longmen Mountains were refugia for F. quadranus during the LGM. Demographic analyses indicated that the Daba lineage experienced population size increase prior to the LGM but the Main Qinling and the Micang–Western Qinling lineages expanded in population size and range after the LGM, and the other lineages almost have stable population size or slight slow population size decline.
Conclusions/Significance
The Qinling–Daba Mountains hosted three refugia for F. quadranus during the LGM. Populations that originated in the Daba Mountains colonized the Main Qinling Mountains after the LGM. Recent sharp expansion of the Micang–Western Qinling and Main Qinling lineages probably contribute to their present-day secondary contact.
Citation: Wang B, Jiang J, Xie F, Li C (2012) Postglacial Colonization of the Qinling Mountains: Phylogeography of the Swelled Vent Frog (Feirana quadranus). PLoS ONE 7(7): e41579. https://doi.org/10.1371/journal.pone.0041579
Editor: Robert DeSalle, American Museum of Natural History, United States of America
Received: March 9, 2012; Accepted: June 22, 2012; Published: July 25, 2012
Copyright: © 2012 Wang et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Funding: This work was supported by the National Natural Sciences Foundation of China (NSFC-31071906, 30730029: http://www.nsfc.gov.cn/Portal0/default152.htm), the Key Laboratory of Mountain Ecological Restoration and Bioresource Utilization in CIB of CAS and Ecological Restoration and Biodiversity Conservation Key Laboratory of Sichuan Province, and the field front project of the knowledge innovation program of CAS ‘Biodiversity conservation and sustainable use of characteristic biological resources’ (Y1C2021200, Y1B3021, KSCX2-EW-J-22:http://www.bio.cas.cn/). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
Climatic changes associated with Pleistocene glacial cycles are believed to have caused montane species to shift, expand, or contract along latitudinal or elevational gradients [1]–[3]. Intraspecific populations adapted to different mountains may have experienced isolations and connectivity as well as range contractions and expansions following climate fluctuations; thus, prolonged isolation may led to substantial genetic heterogeneity among populations and long-term connectivity may contributed to dispersal or ongoing gene flow among mountains; and accompanying with these processes, the population size would be predicted to fluctuate [4]–[6]. The level of divergence together with ecological and biological factors usually affects the formation of independent lineages and whether populations occurring from different regions mix into a same gene pool when they have secondary contact [6]–[9].
Nevertheless, the diversification patterns influenced by Pleistocene climate cycles may have varied from place to place [7]–[9]. Many Europe and North American organisms retreated to the southern regions during the ice ages and recolonized the original north range during the interglacial periods [7]–[12]. In contrast, East Asia was probably characterized by a mosaic of mountains and likely experienced a relatively mild Pleistocene climate [13]–[15], a combination that probably provided stable habitats [16] and glacial refugia for species [17] throughout their entire ranges rather than species being confined only to southern latitudes within their present-day ranges [2], [3], [6]. In consistent to this speculation, some bird species probably persisted in its present range including the Qinling–Daba Mountains and southern China mountains south of Yangtze River during the glacial periods and then expanded during the warming time [18], [19], whereas some herpetology species probably hid into multiple refugia including the south Korea and the regions south of Yangtze River in the southern China during the ice ages [20], [21]; but most of them have probably expanded on population size and range during glacial periods before the last glacial maximum (LGM) [18]–[19]. So far, our knowledge on the influence of the Pleistocene climate changes affecting diversification of the East Asian species is limited to few literatures [18]–[21]; thus, further investigations on phylogeography and demography of the montane species in this region are necessary for testing above distinct hypotheses.
The Qinling–Daba Mountains in the central China comprise of series of mountains, such as Qinling, Daba, Micang, Longmen, and Wuling Mountains, that are delimited by hydrological systems (e.g. Jialing River, Han River and Yangtze River). Moreover, this region supports environments that are nearly unique among those for Oriental and Palearctic organisms that occur in temperate and sub-tropical ecosystems and are inhabited by considerable endemic species [17], [22]. Although most mountains were not glaciated during the Pleistocene [23]–[25], this region experienced climatic fluctuations that likely impacted species distributions, diversification and demography [17]. The frog group Feirana, belonging to the family Dicroglossidae, as stream-dwelling, inhabits this region and is highly dependent on mesic environments that restrict the distributions of them [26]–[31]. Therefore, it was easily affected by the Pleistocene climate changes.
The swelled vent frog Feirana quadranus [26] is restricted to montane streams throughout the Qinling–Daba Mountains [27]–[31]. Because of its distribution, phylogeographic study of this species may contribute to understanding the effects of Pleistocene climate changes on diversification in the Qinling–Daba Mountains. For example, in this species, populations on different mountains could derive from one refugium during the LGM, or may be restricted in multiple refugia and experience long–term vicariances. In the latter speculation, it may be predicted that populations on different mountains could probably merge into a common gene pool during the period of secondary contact. However, based on 12S, 16S and ND2 gene sequences from limited number of samples, Wang et al. [28] indicated that F. quadranus was comprised of several divergent lineages restricted to different mountains and guess that these lineages have possibly experienced prolonged vicariances.
Therefore, in this study, we more intensively explore the phylogeography of F. quadranus based on ND2 gene sequences to test above hypotheses concerning the diversification of this species. First, we use phylogenetic methods to test whether populations on different regions are reciprocally monophyletic. Second, we conduct multiple analyses to test if the identified lineages are isolated by barriers with unsuitable environments that restrict gene flow among them. Third, we use coalescent simulations to test for above alternative refugia hypotheses. Finally, we test whether historical population size fluctuations of lineages were related to the Pleistocene climatic fluctuations.
Materials and Methods
Sampling and Laboratory Protocols
In total, we used 224 samples from 34 populations that cover the entire geographic range of Feirana quadranus (Table 1; Figure 1); data for eighteen samples derives from [28] (GenBank Accession No.: GQ225934, 37, 39, 41–50, 52, 53, 55, 56, 58). Voucher specimens were deposited in Chengdu Institute of Biology (CIB), the Chinese Academy of Sciences (CAS).
Details of the populations are given in Table 1, and the phylogenetic lineages are given in Figure 3.
Total genomic DNA was extracted from muscle tissues (preserved in 95% Ethanol) using a standard phenol-chloroform extraction procedure [32]. To evaluate the presence and extent of mitochondrial introgression between Feirana quadranus and closely related species (F. kangxianensis newly published [33], and F. taihangnica), we amplified the nuclear tyrosinase gene sequences for 68 individuals from 28 localities of F. quadranus, one sample of F. kangxianensis and four samples of F. taihangnica from three localities (Table S1) using primers Tyr1G and Tyr1B [34]. PCR amplification and sequencing were performed according to procedures in other studies [34], [35]. We then obtained mitochondrial ND2 gene sequences for all new samples used in this work according to the methods in our previous study [28]. Following other studies [28], [34], [36], [37], we used Nanorana pleskei, Gyndrapaa yunnanensis, F. kangxianensis, and F. taihangnica as outgroups. Data from these taxa were either sequenced in this present study or downloaded from GenBank (HQ324232, EU979962, EU979976, GQ225989, GQ225994–225996 and GQ226000). New sequences were deposited in GenBank (GenBank Accession No.: JX263025–JX263230, JX263231–JX263298 and JX263299–JX263302).
All sequences were aligned using MEGA 5.0 [38] and manually checked through examining both nucleotide and amino acid sequences to confirm nuclear homologs of target genes. Alignments were unambiguous and no indels were found in these two protein-coding regions. Because the present tyrosinase gene dataset is insufficient for adequately revealing intraspecific relationships (see the results), it was only used to identify whether F. quadranus is an independent evolutionary lineage distinct from F. kangxianensis and F. taihangnica. ND2 gene sequences were used for phylogeographic study of F. quadranus as well as mitochondrial DNA (mtDNA) is commonly used to investigate historical evolutionary patterns [39]–[42].
Phylogenetic Analyses
Haplotypes for our samples were identified by DNASP 5.1 [43], and subsequently used to reconstruct the gene genealogy using maximum parsimony (MP), maximum likelihood (ML) and Bayesian inference (BI), as implemented in PAUP 4.0b10 [44], PHYML 3.0 [45] and MRBAYES 3.0b4 [46], respectively. MP analyses were performed used same parameters as in [28]. For ML and BI analyses, the TrN + G model was selected for all datasets by likelihood ratio tests either under the Corrected Akaike Information Criterion (AICc) [47] or under the Bayesian Information Criterion (BIC) implemented in JMODELTEST 0.1.1 [45], [48]. In ML analyses, the values of the shape parameter of the Gamma distribution were set as 0.186 and 2.032 (estimated by JMODELTEST) for the ND2 and tyrosinase datasets, respectively. The base frequency and ratio of transitions/transversions were optimized by the maximum likelihood criterion in PHYML. The default tree search approach using simultaneous Nearest Neighbor Interchange (NNI) method and BioNJ tree as starting tree was used to estimate tree topologies. Non-parametric bootstrapping with heuristic searches of 10000 replicates was used to assess confidences of branches in MP and ML trees [49]–[52]. In BI, according to the estimations of JMODELTEST, the value of gamma shape was set as 0.186 and the base frequencies was specified as (0.2985 0.3068 0.0954 0.2993); then we initiated two dependent runs each with four simultaneous Monte Carlo Markov chains (MCMC) for 50 million generations with sampling every 1000 generations and discarded the first 10% of generations as “burn-in”. The convergence of chains was confirmed until average standard deviation of split frequency is below 0.01, and the posterior probabilities (pp) were achieved from the remaining trees [53], [54].
Population Structure Analyses
Haplotype network were constructed using parsimony method in the software TCS 1.21 [55] based on ND2 sequences. For mtDNA analyses, DNASP 5.1 was used to compute the number of haplotypes (H), haplotypic diversity (Hd), and nucleotide diversity (π). To test whether there is significant geographical divisions and population genetic structure within this species, an analysis of molecular variance (AMOVA) was implemented in ARLEQUIN 3.1 [56]. The fixation indices Fct, Fsc and Fst that described variation among groups, variation within groups and the variation within populations, respectively, are estimated, and the most probable geographical subdivision were inferred with the highest value of Fct [56]. The significance of these statistics was determined by 10,000 permutation replicates.
To test the isolation by distance model (IBD) with a significant correlation between geographic and genetic distance (Fst/1-Fst), a Mantel test (10,000 randomizations) was performed in the web-based computer program IBDWS 3.14 [57]. Additionally, to highlight geographical discontinuity corresponding to breaks of gene flow among groups and to test if these boundaries are corresponding to major valleys in the Qinling–Daba Mountains (Figure 1), Monmonier’s maximum difference algorithm [58] was implemented in the software BARRIER 2.2 [58]. In this kind of analyses, geographical coordinates were used for each locality and connected by Delauney triangulation using a pairwise Tamura-Nei 93 genetic distance matrix for all populations; then the boundaries (the areas where differences between pairs of populations are largest) were identified using Monmonier’s maximum difference algorithm [58].
Divergence Time Estimates
Under a ‘relaxed molecular clock’ assumption [59], BEAST 1.6.1 [60] was used to estimate the divergence dates among lineages. Under the TrN + G model (estimated by JMODELTEST), genealogy was reconstructed with an uncorrelated lognormal tree prior with an exponential growth prior and a lognormal mean rate (0.65% change per lineage per million years) of ND2 gene which was inferred among previous studies of amphibians and reptiles [61]–[67]. For these analyses, two runs were conducted for 50 million generations with sampling every 1000 generations and in each run, 10% of the initial samples were discarded as burn-in. We checked the convergence of the chains through visual inspection of plotted posterior estimates and by determining the effective sample size (ESS) for each parameter sampled from the MCMC analysis using the program TRACER 1.5 [68]. The final tree including divergence estimates and their 95% highest posterior densities (HPD) were computed in TREEANNOTATOR 1.4.5.
Historical Biogeography
To infer the distributional area of the most recent common ancestor (MRCA) of lineages, we used a ML method conducted in MESQUITE 2.75 [69]. In this analysis, six mountainous regions were regarded as separate character states and each Feirana quadranus individual was coded to one of these regions. The ancestral area was estimated under a ML model on the ML tree. Moreover, the Markov k-state 1 (Mk1) model with equal likelihood migration rates among mountains was selected because there is no prior information for dispersal rates. For each lineage, the area with the highest proportional likelihoods was suggested as the ancestral area.
Statistical phylogeography using coalescent simulations was performed to test the fit of the observed gene trees to phylogeographic hypotheses (see introduction) [8], [70]–[71]. In simulations, the grouping arrangement of populations was based on their geographical origins (Figure 1). Four hypotheses concerning refugia during the last glacial maximum (LGM) were tested for Feirana quadranus [72]. The first is the single-refugium model hypothesizing that all current populations were derived from a single refugium at the end of the LGM (Figure 2A). The second and the third hypotheses predicted that during the LGM there were either two or three refugia for F. quadranus, respectively (Figure 2B, 2C). The final is the multiple-refugia model that posits geographical groups experienced long-term vicariance and were isolated in multiple refugia (Figure 2D). Note that the long-term divergence among groups occurring before the LGM was estimated in this study using BEAST (see the results).
Four hypotheses concerning the refugia during the last glacial maximum (LGM) were tested. a) single-refugium hypothesis, T = 50 ka (LGM is called the Dali glaciation in China); b) two-refugia hypothesis, T = 4.6 Ma, T1 = 50 ka; c) three-refugia hypothesis, T = 4.6 Ma, T1 = 3.3 Ma, T2 = 50 ka; d) multiple-refugia hypothesis, T = 4.6 Ma, T1 = 3.3 Ma, T2 = 2.1 Ma, T3 = 1.8 Ma and T4 = 0.5 Ma. The detail interpretation for these models is given in the text. Branch lengths are time in generations based on a 2.5-year generation time in F. quadranus. Branch widths (effective female population size, Nef) are scaled for each group based on the proportion of the total Nef that each group comprised. Groups A–E occurred from six mountainous regions were given in Figure 1.
Two sets of coalescent simulations were conducted to test hypotheses using the software MESQUITE. In the first set, we simulated 1000 coalescent genealogies constrained within different hypothesized population trees and recorded the distribution of S [73], then through comparing the S value of the empirical ML genealogy with the S values of the simulated genealogies, we tested whether observed genealogies was consistent with the given hypotheses. In the second method, one hundred DNA sequence matrices were simulated under the optimal evolutionary model selected by JMODELTEST for observed dataset, then we reconstructed trees based on these matrices in PAUP and recorded the distribution of S values which were compared with S value of ML genealogy to test model fit. In all simulations, coalescent time (generations) was converted from absolute time (years) divided by 2.5 years as a generation time for Feirana taihangnica [29], [30].
Female effective population size (Nef) of each geographical group for coalescent simulations was converted from Theta using the equation θ = Nefμ with μ = 1.625E-8 (0.65E-8 * 2.5). The θ-values were estimated three times using ML method to ensure convergence in the program MIGRATE-N 3.3 [74] under the following parameters: 10 small chains for 200,000 steps and three long chains for 20 million generations with no heating and chains were sampled every 100 generations following a burn-in of 2 million generations. Total Nef were sum of the Nef for all groups and the proportion of total Nef that each group comprised were scaled the branch width of hypothesized trees (Figure 2) [75]–[77]. The 95% CI’s of Nef calculated from 95% CI’s of θ-values [78] were used as model parameters in simulations.
Historical Demography
Deviations in Tajima’s D [79] and Fu’ Fs [80] implemented in ARLEQUIN were used to test a recent population expansion or bottleneck [79], [80]. The coalescent-based method Bayesian Skyline Plot (BSP) [81], as implemented in the software BEAST, was used to investigate population fluctuations of each lineage over time. The HKY model was selected by JMODELTEST for each lineage. The starting tree was randomly generated. Based on the number of samples, 3–10 grouped coalescent intervals (m) were prior chosen for different lineages; and the piecewise-constant model was selected as a prior skyline model. We conducted two MCMC runs for 20 million iterations, sampling genealogy and population size parameters every 1000 iterations and discarding the first 10% as burn-in. Additionally, we used the mean rate 0.65%/Ma to scale the time axis on BSP and uncorrelated lognormal model to account for rate variation among lineages. Convergence of the chains and ESS was also checked using TRACER.
Results
Phylogenetic Relationships
In the tyrosinase gene dataset (∼521 bps) including outgroup and ingroup sequences, 37 variable sites and 18 polymorphic sites were observed, but in ingroup sequences, we observed only 15 variable sites and 13 polymorphic sites. Twenty-five haplotypes were recognized for tyrosinase gene sequences of Feirana quadranus. Based on tyrosinase gene sequences, F. quadranus was revealed as monophyletic from closely related species (F. kangxianensis and F. taihangnica), but within F. quadranus, the relationships among haplotypes were reticular and not resolved.
An alignment of complete ND2 gene sequences (∼1035 bps) was obtained for all samples and 72 haplotypes were recognized for all Feirana individuals (Table 1). MP analysis produced 1423 maximum parsimonious trees (403 steps; Consistency Index, CI = 0.8488; Retention Index, RI = 0.8036). The log-likelihood score of ML tree was -3097.60666. MP, ML and BI analyses revealed almost consistent topologies except for few tip branches with poor supports. All haplotypes of F. quadranus were clustered into a clade apart from other species (all supports = 100%; Figure 3). Feirana quadranus is composed of three clades including six lineages (Figure 3). The first diverging clade, the Southern Clade, was only comprised of population WL from Wuling Mountains; however, its monophyly receives low supports (<50%; Figures 1, 3). The sister relationship of the Eastern and Western Clades was strongly supported (MP/ML/BI: 82/84/0.95). In the Eastern Clade, the basal lineage (Easternmost Daba lineage) includes only the population CY from the easternmost Daba Mountains, and the Daba lineage included populations from Daba Mountains, whereas the Main Qinling lineage was comprised of individuals from both the Main Qinling and Western Qinling Mountains (Figures 1, 3). Two lineages were found in the Western Clade. One lineage (Longmen lineage) included populations from the Longmen Mountains, and another lineage (Micang–Western Qinling lineage) included most populations from the Micang and Western Qinling Mountains and one population from Main Qinling Mountains (Figures 1, 3). Noticeably, each of populations FX, LD and HX meantime possessed haplotypes of Main Qinling lineage and Micang–Western Qinling lineage (Figure 1).
The bootstrap values calculated by Maximum Parsimony and Maximum Likelihood analyses and the Bayesian posterior probabilities from Bayesian analyses are presented above the main branches; the diamonds on branches mean that all support values are greater than 95%. The estimated divergence time (Ma, mean and 95% CI’s) for major nodes is below the branches. Ten major nodes (1–10) are indicated on the tree. Haplotypes (H1–H72) were given in Table 1.
Population Genetic Structure
TCS analysis based on ND2 gene sequences yielded four unconnected haplotype groups with 95% confidence intervals (connect limit = 14 steps; Figure 4). The Southern Clade connects with the Eastern Clade in 21 steps, and the Eastern Clade connects with the Western Clade in 19 steps. There is no common haplotype among three clades. In the Southern Clade, H1 is connected with H2 and H3 in 15 steps. In the Eastern Clade, the Easternmost Daba lineage connects only with Daba lineage in 11 steps, and the Daba lineage connects with Main Qinling lineage in five steps, but there is no common haplotype or population among these lineages. In the Western Clade, it is need five steps to connect the Micang–Western Qinling lineage with the Longmen lineage. Noticeably, haplotypes H51 and H70 meantime possess individuals from the Western Qinling Mountains and the Main Qinling Mountains and H53 occurring from the Western Qinling Mountains was nested into the Main Qinling lineage.
Haplotypes (H1–H72) correspond to Table 1. The relative sizes of the circles in the network are proportional to the haplotype frequencies (n), and the black dots represent missing haplotypes. The phylogenetic clades and lineages were given in Figure 3.
Substantial Hd was found for all lineages (0.66–0.95), but low π was found for all lineages (0.0016–0.0059) except the Wuling lineage (0.0097; Table 2). We tested the partitioning of genetic variation across four hypothesized sampling groupings according to geographical features and lineage relationships (Table 3; Figures 1, 3). Using AMOVA, significant genetic structure (P<0.01) was found across all hierarchical levels when all group allocations were tested (Table 3). These results rejected the null hypothesis that proposes no genetic structure existing within this species. Highest Fct (0.8117; P<0.001) was achieved when using grouping arrangement based on six phylogenetic lineages, and lowest Fct (0.5243; P<0.001) was showed when grouped all lineages into two large clades. Relatively, when all six lineages were grouped into three clades and five phylogenetic groups, high Fct (0.7182 and 0.7644; P<0.001) were also estimated.
The Mantel test resulted in an uncorrelated relationship (r2 = −0.0784; P>0.05) between genetic distance (Fst/1-Fst) and geological distance. The barrier prediction analysis highlighted several putative barriers to gene flow that separate clades or lineages of Feirana quadranus (Figures 1, 3). The first predicted barrier corresponds to the Qing River that separates the Southern Clade from the Eastern Clade; the second is a barrier basically corresponding to the Qu River–upper Han River–part of the upper Jialing River that separates the Eastern Clade from the Western Clade; the third barrier basically corresponds to the Yangtze River separating the Easternmost Daba lineage from the Daba lineage; the fourth barrier is perfectly in line with the Bailong River separating the Longmen lineage from the Micang–Western Qinling lineage; and the fifth barrier is perfectly in line with the Han River which separates the Main Qinling lineage from the Daba lineage (Figure 1).
Divergence Time, Ancestral Area and Hypotheses Testing
The divergences among lineages within Feirana quadranus were estimated as having occurred from Pliocene to the mid-Pleistocene (Figure 3). The results of the ancestral area reconstruction and divergence dating indicated that F. quadranus probably originated in the Wuling Mountains (Proportional likelihoods is 0.57) during the Miocene (Figure 3; Table 4). The divergence of the Southern Clade from other two major clades possibly occurred about 4.6 Ma (95% CI’s: 0.6–5.9 Ma), and that within this clade maybe occurred in the early Pleistocene (2.2 Ma, 95% CI’s: 0.2–2.4 Ma). The most recent common ancestor (MRCA) of Eastern and Western Clades was most probably in the Daba Mountains (Proportional likelihoods is 0.45), and these lineages approximately diverged about 3.3 Ma (95% CI’s: 0.5–4.5 Ma). Eastern Clade probably originated in the Daba Mountains (Proportional likelihoods is 0.96). In the early Pleistocene (about 2.1 Ma, 95% CI’s: 0.3–2.6 Ma), the Easternmost Daba lineage was possibly separated from the Daba and Main Qinling lineages whose MRCA maybe occurred in the Daba Mountains (Proportional likelihoods is 0.98). The MRCA of the Main Qinling lineage was in the Western Qinling Mountains (Proportional likelihood is 0.99) and maybe diverged from the Daba lineage during the Pleistocene (1.6 Ma, 95% CI’s: 0.2–1.9 Ma). The Western Clade probably originated in its present-day range (Proportional likelihood is 0.37 for each of Micang Mountains and Western Qinling Mountains), and the split between its two lineages (Micang–Western Qinling lineage and Longmen lineage) approximately occurred during the early Pleistocene around 1.8 Ma (95% CI’s: 0.2–1.8 Ma).
Total Nef = 2,043,076 (95% CI: 990,000–5,360,000) was achieved from an overall θ of 0.0332 (95% CI: 0.0162–0.0871) estimated for our ND2 data. The coalescent times of the current trees were simulated under the ‘without migration’ model. The Slatkin and Maddison’s S for our ML genealogy was equal to 21. From both sets of coalescent simulations, the 95% CI’s of S-values of simulated trees (S >23) for the single-refugium model and two-refugia model did not include 21 (S-value of ML genealogy), indicating that these models were rejected (P<0.001), while the multiple-refugia model was also rejected (95% CI’s of S-values <6; P<0.001). However, the three-refugia model could not be rejected (95% CI’s of S-values were 19–28 including 21; P>0.05).
Demographic History
Tajima’s D and Fu’s Fs statistics were positive for the Wuling and Longmen lineages (Table 3) suggesting either constant or contracting population size. The remaining lineages all possess negative Tajima’s D and Fu’s Fs values (Table 3) probably indicating that these lineages have undergone recent expansion. However, Fu’s Fs values were significant only for the Daba and Main Qinling lineages and Tajima’s D values were significant only for the Main Qinling lineage (Table 3). Noticeably, the power of Fu’s Fs was substantially affected by the number of segregating sites although its power was higher than mismatch distribution and Tajima’s D [82]. Therefore, Coalescent methods such as BSP taking genealogy into account may provide a better estimate of demographic history.
The effective sample size (ESS) for each parameter of the BSP was greater than 200, suggesting that the 20 million generations were sufficient to determine the demographic history for each examined lineage. The Wuling and Longmen lineages have maintained slightly slow population decline over time. The Easternmost Daba lineage and the Daba lineage have each experienced slow population growth after about 50,000 and 100,000 years ago, respectively; and at about 20,000 years ago, the Main Qingling and Micang–Western Qinling lineages began sharp expansion on population size. Noticeably, in the Holocene, the Wuling and Longmen lineages probably accelerated decreases in population size, while the Easternmost Daba lineage, Main Qingling lineage, and Micang–Western Qinling lineages maintained stable population sizes and the population size of the Daba lineage decreased. Through comparing recent population size fluctuations among lineages, we found that the Main Qinling lineage displayed the largest increase on population size (about 70-fold) whereas the Micang–Western Qinling lineage experienced about 20-fold population size increase over about past 17,000 years. The population sizes of other lineages changed within 10-fold.
Discussion
Genetic breaks existing among intraspecific populations restricted to different mountains often correspond to low-elevation areas with warmer and drier conditions that form barriers restricting gene flow among populations [7], [12], [72], [76], [77], [83]. Feirana quadranus is restricted to montane streams in six major mountains separated by deep valleys or connected only by narrow passes in the Qinling–Daba Mountains (Figure 1). Accordingly, six divergent lineages structured across these mountains were revealed for F. quadranus (Figures 1, 3 and 4); meanwhile, hierarchical AMOVA (Table 3) and IBD test also indicated restricted gene flow among these lineages. Barrier analyses further highlighted the genetic interruptions between pairs of lineages that occur in the vicinity of the major valleys, such as those of the Qing River, Qu River, Bailong River, Jialing Rive and Han River (Figure 1). Although niche models were not tested for F. quadranus, the existence of putative low-elevation barriers with relative warmer and drier conditions suggest that niche conservatism probably plays an important role in driving diversification within this montane species [84], [85]. Exceptionally, some lineages presently span a short distance across some barriers (e.g., the Main Qinling and Micang–Western Qinling lineages span a short distance across the upper Jialing River and the Daba lineage spans a short distance across the Yangtze River; Figure 1). This supports the idea that these lineages has expanded and contracted their range in the past. The presence of haplotypes that span across barriers (e.g. H51 and H70; Figure 4) is probably attributable to recent migration rather than incomplete lineage sorting [76], [86]. This postulation is also supported by the reality that these haplotypes located in the interior position of phylogenetic tree (Figure 3). If so, it is understandable that the separated effects of some larger rivers (Han River, Jialing River and Yangtze River) are look like lower than that of those big tributaries (Bailong River, Qu River and Qing River). Considerable spatial genetic differentiations among lineages reflect the prolonged historical isolations among them even over Pliocene indicated by coalescent simulations. Thus, although there has been still no direct geological evidence for this phenomenon, we could guess that big tributaries probably have possessed longer geological history than the upriver of those stem river, indicating that these tributaries isolated populations longer than the upper stem rivers.
During Pliocene and Pleistocene, the rapid uplifts of the Tibetan Plateau are hypothesized to have induced environment shifts [87] that may have promoted diversification in the southeastern Asia [88]–[90]. The rapid uplifts of the Tibetan Plateau are believed to correspond to three phases (around 3.6 Ma, 2.6 Ma and 1.7 Ma, respectively) that are accompanied by the formation of the Asian monsoon, the beginning of Chinese loess, and appearance of the Yellow River [87], [91], although the time of origin of the current Asian monsoon is still debated [92], [93]. The geographical feature of Qinling–Daba Mountains at the eastern portion of the Tibetan Plateau was affected by these drastically movements [87], [94]. Divergence time estimations and coalescent simulations suggested that the divergence among three major clades within Feirana quadranus occurred approximately during this period (Figures 2 and 3). Consequently, we suggested that the basal diversification of this species were probably affected by the significant tectonic events associated with the rapid uplift of the Tibetan Plateau during the Pliocene–early Pleistocene.
Historical allopatry and multiple refugia during the LGM have been suggested by previous studies of East Asian animals [18]–[21]. Through coalescent simulations for Feirana quadranus, we found three refugia located in the Wuling Mountains, Daba Mountains, and the Longmen–Micang–Western Qinling Mountains, respectively (Table 4; Figure 2C). The first two are consistently recognized as refugia for frogs [20], reptiles [21], and birds [19]. But the third region was firstly exactly recognized as a refugium for vertebrates using an effective method (coalescent simulations), although few studies indicated this point with no direct result as evidence [18], [19]; whereas, this is in accordance with expectations. As mentioned above, these regions possessed suitable environments for many temperate and subtropical species [15], [22]. During glacial periods, local climate changes in this region were probably ameliorated by orographic precipitation caused by the windward slopes of the Tibetan Plateau and Qinling Mountains though it is inland with weaker summer Asia monsoons. The moist condition in these regions during the last glacial period is also supported by palaeoclimatic evidence [14], [95]. It is conceivable that in these refugia F. quadranus survived in the lower elevations during the LGM, with subsequent dispersals (Figure 5). Of course, additional cases with congruent phylogeographic patterns will further reveal the legacies of glacial cycles for organisms in this region. Additionally, in contrast to phylogeographical studies on birds [18], [19] but consistent with those on the black-spotted frog [20] and short-tailed pit viper [21], we did not recover the Main Qinling Mountains as a refugium for the swelled vent frog. Instead, based on coalescent simulations and ancestral area reconstructions, we propose that populations presently holding this region originated from the Daba Mountains and colonized the Main Qinling Mountains and even some part of the Western Qinling Mountains after the LGM (Table 4). Colder conditions in the Main Qinling Mountains, including ice sheets in some mountains (e.g. Taibai Mountain) during the LGM [96], [97], possibly excluded this poikilotherm from this region, which originated in more southern areas. Moreover, the competition from its sibling species F. taihangnica which is native to the Qinling Mountains [28] had also probably contributed to the absence of this species in the Main Qinling Mountains during the ice ages.
The central solid line represents the median value for the log of the population size (Nef * τ) and the area between two thinner lines represent the 95% higher posterior density. The thicker dashed line represents the median of estimation time of the most recent common ancestor (MRCA), and the thinner is the lower limit of the 95% confidence interval. The combined figure shows the population size fluctuation after the lower limit of 95% CI’s of estimation time of the most recent common ancestor. Phylogenetic lineages correspond to Figure 3.
The dispersal model also illustrates the demographic history and historical vicariance of this frog species. As one example, the secondary contact and overlapping ranges of haplotypes occurred between the Main Qinling lineage and the Micang–Western Qinling lineage is supported (Figures 3, 4). Two haplotypes (H51 and H53) representing individuals from four populations in the Western Qinling Mountains were resolved as parts of the Main Qinling lineage, as well as some individuals of haplotype H70 now occurs from the Main Qinling Moutains was nested into the Micang–Western Qinling lineage (Figure 4). The expansion and secondary contact presumably occurred after the LGM because of the relative interior position of these lineages and based on coalescent simulations. This hypothesis is supported by BSP analyses which reveal that the Main Qinling and Micang–Western Qinling lineages experienced sharp expansions during the present interglacial period in China (Figure 5) [98]. As another example, ancestral area reconstructions indicated that the ancestor of the Eastern Clade was in the Daba Mountains (Table 4), inferring that populations in the Main Qinling Mountains have dispersed from the Daba Mountains and invaded the Qinling Mountains after the LGM. This idea is also supported by the coalescent simulations and the polyphyletic relationships among lineages within the Eastern Clade (Figure 3). Based on the topography of this region, we suggest the hilly region north of the population ZB as the dispersal conduit for the ancestor of the Main Qinling lineage from the Daba Mountains to the Main Qinling Mountains, which is geographically distinct from the lower-elevation valley Han River (Figure 1). Currently, the absence of this frog in the vicinity of the Han River indicates that gene flow between the Daba and Main Qinling lineages has decreased. On the other hand, the biogeographic events of this species reflect that the Han River probably ran through the hilly region between Daba and Qinling Mountains after the LGM.
According to our mtDNA and nuclear genealogies, we found no hybridization between F. quadranus and its sister group, but we could not determine whether extensive gene flow via male dispersal occurred between mtDNA lineages within F. quadranus because nuclear sequences examined here provided only a small amount of information. If there is hybridization within F. quadranus, it is likely that lineages hybridize within the areas where they come into secondary contact, yet this small amount of hybridization is not enough to obscure historic effects of prolonged isolation among lineages [99], [100], and mitochondrial introgression likely does not extend much beyond contact zones. Still, it is desirable to examine extensive gene flow via male dispersal between lineages [101], [102] through deep studies based on multiple nuclear loci, which may either corroborate or overturn the observed mtDNA phylogenographic patterns of this species.
The repetitive geographic shifts and population size fluctuations in species associated with Pleistocene glacial oscillations are expected to be reflected by increases or decreases in levels of genetic variation and coalescence times [1]–[3]. In Europe and North America, species expanded from southern refugia after the LGM [103], and the genetic diversity of these species is usually geographically structured with declines in genetic diversity towards the north [1]. In this study, genetic diversity did not show a significant north–south trend among populations of Feirana quadranus, suggesting the existing of multiple refugia along the Qinling and Daba Mountains during the LGM. The significantly negative Fu’s Fs suggest recent expansion of the Main Qinling and Daba lineages and the BSPs consistently revealed population growth curves for these lineages, (Table 3; Figure 5). Differently, BSPs showed recent population increase but without significantly negative Fu’s Fs and Tajima’s D for Easternmost Daba and Micang–Western Qinling lineages. The interpretation for this discordance of neutral tests and BSPs could mean that populations underwent recent expansion but then were recently subdivided, subjected to substantial migration, and/or had experienced historical contractions [104]–[106]. Noticeably, the highest local genetic diversity was found in the Southern Clade (Wuling lineage) that includes only one population (Table 3). This could be partially attributed to historical contraction or subdivision evidenced by positive Tajima’s D and Fu’s Fs values [79], [80]. Similarly, the Longmen lineage exhibits relative low genetic diversity and positive Fu’s Fs, and Tajima’s D values (Table 3) that indicate that this lineage probably had experienced bottlenecks [79], [80] which was also supported by the BSP of this lineage (Figure 5). Of course, the limitation of sampling for these lineages probably affected our results; so we take a caution to make these speculations and would include more information such as sequences and genes to test these hypotheses in the future.
Some biological factors of amphibians, such as higher sensitivity to climatic changes [107], [108] and poor dispersal capability, probably restricted them into different refugia and experienced contractions during the LGM. The slight slow decline of the Wuling and Longmen lineages supports this speculation although BSPs could not deduce more long time coalescent events (Figure 5). The times of expansion among lineages look like different (Figure 5). The Daba and Easternmost Daba lineages probably expanded prior to the LGM (Figure 5) similar to the patterns observed for reptile and bird species in the southeastern China [18], [19], [21]. The interpretation of these previous studies was that the moderate climate after the Marine Isotope Stages (MIS) 16–18 [109] probably provided environments persistence through the glaciations. These studies suggested that the region in the southern part of Daba Mountains might be one refugium for these species [19]–[21], which also presumably supplied conditions for expansion of the Daba and Easternmost Daba lineages prior to the LGM. However, the sharp population size increases detected for the Main Qinling and Micang-Western Qinling lineages occurred after about 20 ka which is roughly consistent with the LGM during MIS2. Conceivably, the MRCA of Main Qinling lineage expanded significantly from the Daba Mountains at the end of the LGM. This speculation is also indicated by the coalescent simulations. As remarked above, during the glacial age, unfavorable environments in the Main Qinling Mountains north of the Daba Mountains and competition from closely related species had possibly restrained populations of F. quadranus from colonizing the Main Qinling Mountains.
Supporting Information
Table S1.
Sampling information and haplotypes of Feirana based on nuclear tyrosinase gene.
https://doi.org/10.1371/journal.pone.0041579.s001
(XLS)
Acknowledgments
We are grateful to Liang Gang, Chen Xiaohong, Hu Junhua, Wang Jie, and Yang Xin for help in collecting samples, Dr. David C. Blackburn (California Academy of Sciences) for polishing English.
Author Contributions
Conceived and designed the experiments: JJ BW. Performed the experiments: BW. Analyzed the data: BW. Contributed reagents/materials/analysis tools: BW FX CL. Wrote the paper: BW JJ FX CL.
References
- 1. Hewitt GM (1996) Some genetic consequences of ice ages, and their role in divergence and speciation. Biol J Linn Soc 58: 247–276.
- 2. Hewitt GM (2000) The genetic legacy of the Quaternary ice ages. Nature 405 (6789): 907–913.
- 3. Hewitt GM (2004) Genetic consequences of climatic oscillations in the Quaternary. Philos T R Soc B 359(1442): 183–195.
- 4. Avise JC, Walker D (1998) Pleistocene phylogeographic effects on avian populations and the speciation process. P Roy Soc B-Biol Sci 265(1395): 457–463.
- 5. Avise JC, Walker D, Johns GC (1998) Speciation durations and Pleistocene effects on vertebrate phylogeography. P Roy Soc B-Biol Sci 265(1407): 1707–1712.
- 6.
Avise JC (2000) Phylogeography: The History and Formation of Species. Cambridge, Harvard University Press.
- 7. Knowles LL (2000) Tests of Pleistocene speciation in montane grasshoppers (genus Melanoplus) from the sky islands of western North America. Evolution 54: 1337–1348.
- 8. Knowles LL (2001) Did the Pleistocene glaciations promote divergence? Tests of explicit refugial models in montane grasshoppers. Mol Ecol 10: 691–701.
- 9. Maddison WP, McMahon M (2000) Divergence and reticulation among montane populations of a jumping spider (Habronattus pugillis Griswold). Syst Biol 49: 400–421.
- 10. Masta SE (2000) Phylogeography of the jumping spider Habronattus pugillis (Araneae: Salticidae): recent vicariance of sky island populations? Evolution 54: 1699–1711.
- 11. Carstens BC, Brunsfeld SJ, Demboski JR, Good JM, Sullivan J (2005) Investigating the evolutionary history of the Pacific Northwest mesic forest ecosystem: hypothesis testing within a comparative phylogeographic framework. Evolution 59: 1639–1652.
- 12. Smith CI, Farrell BD (2005) Phylogeography of the longhorn cactus beetle Moneilema appressum LeConte (Coleoptera: Cerambycidae): was the differentiation of the Madrean sky islands driven by Pleistocene climate changes? Mol Ecol 14: 3049–3065.
- 13. Pinot S, Ramstein G, Harrison SP, Prentice IC, Guiot J, et al. (1999) Tropical paleoclimates at the last glacial maximum: comparison of Paleoclimate Modeling Intercomparison Project (PMIP) simulations and paleodata. Clim Dynam 15: 857–874.
- 14. Ju L, Wang H, Jiang D (2007) Simulation of the Last Glacial Maximum climate over East Asia with a regional climate model nested in a general circulation model. Palaeogeogr Palaeocl 248: 376–390.
- 15.
Zhang RZ (1999) Zoogeography of China. Beijing, Science Press.
- 16. Qian H, Ricklefs RE (2000) Large-scale processes and the Asian bias in temperate plant species diversity. Nature 407: 180–182.
- 17. Zhang RZ (2004) Relict distribution of land vertebrates and Quaternary glaciation in China. Acta Zool Sinica 50: 841–851.
- 18. Song G, Qu YH, Yin ZH, Li SS, Liu NF, et al. (2009) Phylogeography of the Alcippe morrisonia (Aves: Timaliidae): long population history beyond late Pleistocene glaciations. BMC Evol Biol 9: 143.
- 19. Li SH, Yeung CKL, Feinstein J, Han L, Le MH, et al. (2009) Sailing through the Late Pleistocene: unusual historical demography of an East Asian endemic, the Chinese Hwamei (Leucodioptron canorum canorum), during the last glacial period. Mol Ecol 18(4): 622–633.
- 20. Zhang H, Yan J, Zhang GQ, Zhou KY (2008) Phylogeography and demographic history of chinese black-spotted frog Populations (Pelophylax nigromaculata): Evidence for independent refugia expansion and secondary contact. BMC Evol Biol 8: 21.
- 21. Ding L, He SP, Gan XN, Zhao EM (2011) A phylogeographic, demographic and historical analysis of the short-tailed pit viper (Gloydius brevicaudus): Evidence for early divergence and late expansion during the Pleistocene. Mol Ecol 20: 1905–1922.
- 22.
Chen L (2008) The Precise Biogeographic Division of Palaearctic and Oriental Realm in the East of China: Based on Data of Amphibians [D]. Northwest Institute of Plateau Biology, the Chinese Academy of Sciences.
- 23. Shi YF (2002) Characteristics of late Quaternary monsoonal glaciation on the Tibetan Plateau and in East Asia. Quatern Int 97/98: 79–91.
- 24.
Shi YF (2007) The Pleistocene Glaciation of China and Environmental Change. Zhengzhou, Henan Science and Technology Publishing House.
- 25. Li JJ, Shu Q, Zhou SZ, Zhao ZJ, Zhang JM (2004) Review and prospects of Quaternary glaciation research in China. J Glaciol Geocryol 26 (3): 235–243.
- 26. Liu CZ, Hu SC, Yang FH (1960) Amphibians from Wushan, Szechwan. Acta Zool Sin 2: 278–292.
- 27. Wang B, Jiang JP, Chen XH, Xie F, Zheng ZH (2007) Morphometrical study on populations of the genus Feirana (Amphibia, Anura, Ranidae). Acta Zoot Sin 3: 139–145.
- 28. Wang B, Jiang JP, Xie F, Chen XH, Dubois A, et al. (2009) Molecular Phylogeny and Genetic Identification of Populations of Two Species of Feirana Frogs (Amphibia: Anura, Ranidae, Dicroglossinae, Paini) Endemic to China. Zool Sci 26: 500–509.
- 29.
Fei L, Ye CY, Huang YZ, Jiang JP, Xie F (2005) An Illustrated Key to Chinese Amphibians. Sichuan Publishing House of Science and Technology, Chengdu.
- 30.
Fei L, Hu SQ, Ye CY, Huang YZ (2009) Fauna Sinica, Amphibia: Vol 3. Anura [M]. Beijing: Science Press.
- 31.
Fei L, Ye CY, Jiang JP (2010) Colored Atlas of Chinese Amphibians. Sichuan Publishing House of Science and Technology, Chengdu.
- 32.
Hillis DM, Mable BK, Larson A, Davis SK, Zimmer EA (1996) Nucleic acids IV: sequencing and cloning. In “Molecular Systematics”. Ed by D Hillis, C Moritz, B Mable, Sinauer Associates, Sunderland, Massachusetts, USA, 321–378.
- 33. Yang X, Wang B, Hu JH, Jiang JP (2011) A New Species of the Genus Feirana (Amphibia: Anura: Dicroglossidae) from the Western Qinling Mountains of China. Asian Herpetol Res 2(2): 72–86.
- 34. Che J, Hu JS, Zhou WW, Murphy RW, Papenfuss TJ, et al. (2009) Phylogeny of the Asian spiny frog tribe Paini (Family Dicroglossidae) sensu Dubois. Mol Phylogenet Evol 50: 59–73.
- 35. Chiari Y, Vences M, Vieites DR, Rabemananjara F, Bora P, et al. (2004) New evidence for parallel evolution of colour patterns in Malagasy poison frogs (Mantella). Mol Ecol 13: 3763–3774.
- 36. Jiang JP, Dubois A, Ohler A, Tillier A, Chen XH, et al. (2005) Phylogenetic relationships of the tribe Paini (Amphibia, Anura, Ranidae) based on partial sequences of mitochondrial 12 s and 16 s rRNA genes. Zool Sci 22: 353–362.
- 37. Che J, Zhou WW, Hu JS, Yan F, Papenfuss TJ, et al. (2010) Spiny frogs (Paini) illuminate the history of the Himalayan region and Southeast Asia. P Natl Acad Sci USA 107(31): 13765–13770.
- 38. Tamura K, Peterson D, Peterson N, Stecher G, Nei M, et al. (2011) MEGA5: Molecular Evolutionary Genetics Analysis using Maximum Likelihood, Evolutionary Distance, and Maximum Parsimony Methods. Mol Biol Evol 28: 2731–2739.
- 39. Eizirik E, Bonatto SL, Johnson WE, Crawshaw PG, Vié JC, et al. (1998) Phylogeographic patterns and evolution of the mitochondrial DNA control region in two Neotropical cats (Mammalia, Felidae). J Mol Evol 47: 613–624.
- 40. Culling MA, Janko K, Boron A, Vasil’ev VP, Cöté IM, et al. (2006) European colonization by the spined loach (Cobitis taenia) from Ponto-Caspian refugia based on mitochondrial DNA variation. Mol Ecol 15: 173–190.
- 41. Tchaicka L, Eizirik E, De Oliveira TG, Candido JF, Freitas TRO (2007) Phylogeography and population history of the crab-eating fox (Cerdocyon thous). Mol Ecol 16: 819–838.
- 42. Hudson RR, Turell M (2003) Stochasiticity overrules the “three-times rule”: genetic drift, genetic draft, and coalescence times for nuclear loci versus mitochondrial DNA. Evolution 57: 182–90.
- 43.
Librado P, Rozas J (2009) DnaSP v5: A software for comprehensive analysis of DNA polymorphism data. Bioinformatics 25: 1451–1452 |. https://doi.org/10.1093/bioinformatics/btp187
- 44.
Swofford D (2002) PAUP*. Phylogenetic analysis using parimony (* and other methods). 4.0th edition. Sunderland, Massachusetts: Sinauer Associates.
- 45. Guindon S, Gascuel O (2003) A simple, fast, and accurate algorithm to estimate large phylogenies by maximum likelihood. BMC Syst Biol 52: 696–704.
- 46. Ronquist FR, Huelsenbeck JP (2003) MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics 19: 1572–1574.
- 47. Hurvich CM, Tsai CL (1989) Regression and time series model selection in small samples. Biometrika 76: 297–307.
- 48. Posada D (2008) jModelTest: Phylogenetic model averaging. Mol Biol Evol 25(7): 1253–1256.
- 49. Felsenstein J (1985) Confidence limits on phylogenies: an approach using the bootstrap. Evolution 39: 783–791.
- 50. Felsenstein J, Kishino H (1993) Is there something wrong with the bootstrap on phylogeny? A reply to Hillis and Bull. Syst Zool 42: 193–200.
- 51. Hedges SB (1992) The number of replications needed for accurate estimation of the bootstrap P-value in phylogenetic studies. Mol Biol Evol 9: 366–369.
- 52. Huelsenbeck JP, Hillis DM (1993) Success of phylogenetic methods in the four-taxon case. Syst Biol 42: 247–264.
- 53. Leaché AD, Reeder TW (2002) Molecular systematics of the eastern fence lizard (Sceloporus undulatus): a comparison of parsimony, likelihood, and Bayesian approaches. Syst Biol 51: 44–68.
- 54. Parra-Olea G, García-París M, Wake DB (2004) Molecular diversification of salamanders of the tropical American genus Bolitoglossa (Caudata: Plethodontidae) and its evolutionary and biogeographical implications. Biol J Linn Soc 81: 325–346.
- 55. Clement M, Posada D, Crandall KA (2000) TCS: a computer program to estimate gene genealogies. Mol Ecol 9 (10): 1657–1659.
- 56. Excoffier L, Laval G, Schneider S (2005) Arlequin version 3.0: an integrated software package for population genetics data analysis. Evol Bioinform Online 1: 47–50.
- 57. Jensen JL, Bohonak AJ, Kelley ST (2005) Isolation by distance, web service v.3.14. BMC Genetics 6: 13.
- 58. Manni F, Guerard E, Heyer E (2004) Geographic patterns of (Genetic, Morphologic, Linguistic) Variation: How Barriers Can Be Detected by Using Monmonier’s Algorithm. Hum Biol 76(2): 173–190.
- 59. Drummond AJ, Ho SYW, Phillips MJ, Rambaut A (2006) Relaxed phylogenetics and dating with confidence. Plos Biol 4: 699–710.
- 60. Drummond AJ, Rambaut A (2007) BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol Biol 7: 214.
- 61.
Bermingham E, McCafferty SS, Martin AP (1997) Fish biogeography and molecular clocks: Perspectives from the Panamanian Isthmus. Pp, 113–118 in Kocher TD and Stepien CA, (eds.), Molecular Systematics of Fishes. San Diego: Academic Press.
- 62. Weisrock DW, Macey JR, Ugurtas IH, Larson A, Papenfuss TJ (2001) 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, Mol Phylogenet Evol 18: 434–448.
- 63. Macey JR, Schulte JAII, Larson A, Fang Z, Wang Y, et al. (1998a) Phylogenetic relationships of toads in the Bufo bufo species group from the eastern escarpment of the Tibetan Plateau: A case of vicariance and dispersal. Mol Phylogenet Evol 9: 80–87.
- 64. Macey JR, Schulte JAII, Ananjeva NB, Larson A, Rastegar-Pouyani N, et al. (1998b) Phylogenetic Relationships among Agamid Lizards of the Laudakia caucasia Species Group: Testing Hypotheses of Biogeographic Fragmentation and an Area Cladogram for the Iranian Plateau. Mol Phylogenet Evol 10: 118–131.
- 65. Macey JR, Wang Y, Ananjeva NB, Larson A, Papenfuss TJ (1999) Vicariant patterns of fragmentation among gekkonid lizards of the genus Teratoscincus produced by the Indian Collision: A molecular phylogenetic perspective and an area cladogram for Central Asia. Mol Phylogenet Evol 12: 320–332.
- 66. Macey JR, Strasburg J, Brisson J, Vredenburg V, Jennings M, et al. (2001) Molecular phylogenetics of western North American frogs of the Rana boylii species group. Mol Phylogenet Evol 19: 131–143.
- 67. Moritz C, Dowling TE, Brown WM (1987) Evolution of animal mitochondrial DNA: Relevance for population biology and systematics. Annu Rev Ecol Evol S 18: 269–292.
- 68.
Rambaut A, Drummond AJ (2007) Accessed 2011 Aug 18. Tracer v1.4, Available: http://beast.bio.ed.ac.uk/Tracer.
- 69.
Maddison WP, Maddison DR (2008) Mesquite: a modular system for evolutionary analysis. Version 2.5. Accessed 2011 September 30.
- 70. Knowles LL, Maddison WP (2002) Statistical phylogeography. Mol Ecol 11: 2623–2635.
- 71. Richards CL, Carstens BC, Knowles LL (2007) Distribution modelling and statistical phylogeography: an integrative framework for generating and testing alternative biogeographical hypotheses. J Biogeogr 34: 1833–1845.
- 72. Carstens BC, Knowles LL (2007) Shifting distributions and speciation: Species divergence during rapid climate change. Mol Ecol 16: 619–627.
- 73. Slatkin M, Maddison WP (1989) A cladistic measure of gene flow inferred from phylogenies of alleles. Genetics 123: 603–613.
- 74.
Beerli P (2008) Accessed 2011 September 1. MIGRATE Version 2.4.1, Available: http://popgen.csit.fsu.edu/Migrate-n.html.
- 75. Carstens BC, Sullivan J, Davalos LM, Larsen PA, Pedersen SC (2004) Exploring population genetic structure in three species of Lesser Antillean bats. Mol Ecol 13: 2557–2566.
- 76. Shepard DB, Burbrink FT (2008) Lineage diversification and historical demography of a sky island salamander, Plethodon ouachitae, from the Interior Highlands. Mol Ecol 17: 5315–5335.
- 77. Shepard DB, Burbrink FT (2009) Phylogeographic and demographic effects of Pleistocene climatic fluctuations in a montane salamander, Plethodon fourchensis. Mol Ecol 18: 2243–2262.
- 78. Edwards SV, Beerli P (2000) Gene divergence, population divergence, and the variance in coalescence time in phylogeographic studies. Evolution 54: 1839–1854.
- 79. Tajima F (1989) Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics 123(3): 585–595.
- 80. Fu YX (1997) Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics 147: 915–925.
- 81. Drummond AJ, Rambaut A, Shapiro B, Pybus OG (2005) Bayesian coalescent inference of past population dynamics from molecular sequences. Mol Biol Evol 22(5): 1185–1192.
- 82. Ramos-Osins SE, Rozas J (2002) Statistical properties of new neutrality tests against population growth. Mol Biol Evol 19: 2092–2100.
- 83. DeChaine EG, Martin AP (2006) Using coalescent simulations to test the impact of Quaternary climate cycles on divergence in an alpine plant–insect association. Evolution 60: 1004–1013.
- 84. Wiens JJ, Graham CH (2005) Niche conservatism: integrating evolution, ecology, and conservation biology. Annu Rev Ecol Evol S 36: 519–539.
- 85. Kozak KH, Wiens JJ (2006) Does niche conservatism promote speciation? A case study in North American salamanders. Evolution 60: 2604–2621.
- 86. Peters JL, Zhuravlev Y, Fefelov I, Logie A, Omland KE (2007) Nuclear loci and coalescent methods support ancient hybridization as cause of mitochondrial paraphyly between gadwall and falcated duck (Anas spp.). Evolution 61(8): 1992–2006.
- 87. Li JJ, Fang XM, Pan BT, Zhao ZJ, Song YG (2001) Late Cenozoic intensive uplift of Qinghai-Xizang Plateau and its impacts on environments in surrounding area. Quaternary Sci 21(5): 381–391.
- 88.
Xu R (1981) On palaeobotanical evidence for continental drift and the Himalayan uplift [A]. In: Academia Sinica Tibet Expedition Team, editor. Beijing: Science Press 8–18.
- 89. Guo X, He S, Zhang Y (2005) Phylogeny and biogeography of Chinese sisorid catfishes re-examined using mitochondrial cytochrome b and 16S rRNA gene sequences. Mol Phylogenet Evol 35: 344–362.
- 90. Guo P, Liu Q, Li C, Cheng X, Jiang K, et al. (2011) Molecular phylogeography of Jerdon’s pitviper (Protobothrops jerdonii): importance of the uplift of the Tibetan plateau. J Biogeogr 38: 2326–2336.
- 91. Zhong DL, Ding L (1996) The discuss mechanism and uplift process of Qinghai-Tibet plateau [J]. Sci China Earth Sci 26(4): 289–295.
- 92. An ZS, Kutzbach JE, Prell WL, Porter SC (2001) Evolution of Asian monsoons and phased uplift of the Himalaya-Tibetan plateau since Late Miocene times. Nature 411: 62–66.
- 93. Molnar P (2005) Mio-Pliocene growth of the Tibetan Plateau and evolution of East Asian climate. Palaeontol Electron 8(1): l–23.
- 94. Pan BT, Gao HS, Li BY, Li JJ (2004) Step-like landform sand up lift of the Qinghai-Xizang Plateau. Quaternary Sci 24(1): 50–57.
- 95. Knutti R, Flückiger J, Stocker TF, Timmermann A (2004) Strong hemispheric coupling of glacial climate through freshwater discharge and ocean circulation. Nature 430: 851–856.
- 96.
Shi YF, Li JJ, Li BY (1998) Uplift and environmental changes of Qinghai–Tibetan plateau in the Late Cenozoic. Guangdong Science and Technology Press, Guangzhou.
- 97. Shi YF, Li JJ, Li BY, Yao TD, Wang SM, et al. (1999) Uplift of the Qinghai-Xizang (Tibetan) Plateau and East Asia enviromental change during late cenozoic. Acta Geographica Sinica 54: 10–20.
- 98. Cui ZJ, Zhang W (2003) Discussion about the glacier extent and advance/retreat asynchrony during the last glaciation. J Glaciol Geocryol 25: 510–516.
- 99. Weisrock DW, Kozak KH, Larson A (2005) Phylogeographic analysis of mitochondrial gene flow and introgression in the salamander, Plethodon shermani. Mol Ecol 14: 1457–1472.
- 100. Weisrock DW, Larson A (2006) Testing hypotheses of speciation in the Plethodon jordani species complex with allozymes and mitochondrial DNA sequences. Biol J Linn Soc 89: 25–51.
- 101. Jockusch EL, Wake DB (2002) Falling apart and merging: diversification of slender salamanders (Plethodontidae: Batrachoseps) in the American West. Biol J Linn Soc 76: 361–391.
- 102. Keogh JS, Webb JK, Shine R (2007) Spatial genetic analysis and long-term mark-recapture data demonstrate male-biased dispersal in a snake. Biol Letters 3: 33–35.
- 103. Rowe KC, Heske EJ, Brown PW, Paige KN (2004) Surviving the ice: northern refugia and postglacial colonization. P Natl Acad Sci USA 101: 10355–10359.
- 104. Ray N, Currat M, Excoffier L (2003) Intra-deme molecular diversity in spatially expanding populations. Mol Biol Evol 20: 76–86.
- 105. Castoe TA, Spencer CL, Parkinson CL (2007) Phylogeographic structure and historical demography of the western diamondback rattlesnake (Crotalus atrox): A perspective on North American desert biogeography. Mol Phylogenet Evol 42: 193–212.
- 106. Bertorelle G, Slatkin M (1995) The number of segregating sites in expanding human populations, with implications for estimates of demographic parameters. Mol Biol Evol 12: 887–892.
- 107.
Carey C, Bradford DF, Brunner JL, Collins JP, Davidson EW, et al. (2003) Biotic factors in amphibian population declines. In: Linder G, Krest SK, Sparling DW, editors. pp. 153–208. Pensacola: Society of Environmental Toxicology and Chemistry Press.
- 108. Collins JP, Storfer A (2003) Global amphibian declines: Sorting the hypotheses. Divers Distrib 9: 89–98.
- 109. Wu GJ, Pan BT, Guan QY, Gao HS (2002) The maximum glaciation and desert expansion in China during MIS16. J Glaciol Geocryol 24: 544–549.