Population Structure, Historical Biogeography and Demographic History of the Alpine Toad Scutiger ningshanensis in the Tsinling Mountains of Central China

Population genetic structure, historical biogeography and historical demography of the alpine toad Scutiger ningshanensis were studied using the combined data mtDNA cytochrome b (cyt b) and the mtDNA cytochrome c oxidase subunit I (COI) as the molecular markers. This species has high genetic variation. There was a significant genetic differentiation among most populations. Three lineages were detected. The phylogenetic relationship analyses and the SAMOVA (spatial analysis of molecular variance) results showed significant phylogeographic structure. 82.15% genetic variation occurred among populations whereas differentiation within populations only contributed 17.85% to the total. Mantel test results showed a significant correlation between the pairwise calculated genetic distance and pairwise calculated geographical distance of the populations (regression coefficient  = 0.001286, correlation coefficient  = 0.77051, p (rrand≥robs)  = 0.0185<0.05), indicating the existence of isolation-by-distance pattern of genetic divergence for cyt b + COI sequence, which suggests that the distribution of genetic variation is due to geographical separation rather than natural selection. The population expansion or contraction and genetic differentiation between populations or lineages could be explained by topography and the repetitive uplifts of the Tsinling Mountains and the climatic cycles during the late Pliocene and Pleistocene. S. ningshanensis experienced a rapid population expansion about 40,000 years before present. The current decline in population size was probably caused by anthropogenic disturbance. Current populations of S. ningshanensis are from different refugia though the location of these refugia could not be determined in our study. Topography, climatic changes and repetitive population expansion/contraction together led to the high level of genetic variation in S. ningshanensis. A total of three management units (MUs) was determined, which must be considered when conservation policy is made in the future.


Introduction
Population genetic structure refers to the geographical pattern of genetic diversity within or among populations. It could be influenced by gene flow, genetic drift, selection, mutation and recombination. Gene flow is caused by the movement of individuals from one population to another [1]. Estimation of the gene flow level allows conservation biologists to understand the relationships between populations and assess levels of genetic variation in order to evaluate the relative levels of conservation concern hierarchically across populations in a species. Genetic drift is the change in the frequency of a gene variant in a population due to random sampling [2]. http://en.wikipedia.org/ wiki/Genetic_drift -cite_note-Masel_2011-1#cite_note-Ma-sel_2011-1Genetic drift may lead to disappearance of gene variants and thereby reduce genetic diversity.
Phylogeography connects historical processes in evolution with spatial distributions [4]. Analysis of mitochondrial data promoted the empirical development of phylogeography [3]. The statistical phylogeography is one of the widely used approaches in phylogeography, which takes into account the stochasticity of genetic processes into demographic inference based on coalescent models for parameter estimation [4,5].
The Tsinling Mountains are located in the central part of China, stretching from west to east (Fig. 1). These Mountains are boundary between Oriental realm and Palaearctic realm according to the zoogeographical regions of China [6], and also the watershed for Yangtze River and Yellow River catchment areas, for climate, flora and fauna in China [7][8][9]. The Oriental and Palaearctic species congregate here forming a specific biotic province and containing rich animal and plant resources [10,11]. Like other regions in the northern hemisphere, the Tsinling Mountains experienced several glacial-interglacial cycles during Pleistocene [12][13][14][15][16][17][18]. The climate associated with Pleistocene glacial cycles in East Asia was likely mild and characterized by a mosaic of mountains [19,20]. The past climatic events, such as the Quaternary glaciation, are believed to have played an important role in forming the geographical pattern of the montane species and could leave the vestiges in geographical distribution of genetic diversity of population [21][22][23][24][25][26]. The founder effect during the postglacial population recovery causes a reduction in population genetic diversity [27,28], and the subsequent rapid population expansion [29] may erase the previous geographical differences of the genetic diversity.
The alpine toad Scutiger ningshanensis was described from the western part of the Tsinling Mountains [30] (Fig. 1). Four years later, the second specimen of this species was collected from the same locality [31]. Since then, other specimen of this species was not collected until 2009 when some specimens were collected from several localities in the eastern part of the Tsinling Mountains [32]. Other than the reports on collection of additional specimens, only the biological characteristics of tadpoles of this species was studied [33]. The habitat of this species was roughly divided into two parts: the western part and the eastern part. Is this geographical pattern caused by habitat fragmentation or by populations from different glacial refugia? Does the isolation by distance between the local populations result in occurrence of any speciation events? The aims of the present study were to explore the population genetic structure, historical biogeography and the historical demography of S. ningshanensis.

Ethics statement
This study was approved by the Institutional Animal Care and Use Committee (IACUC) of Shaanxi Normal University. Only the clipped toes or tail tips of tadpoles were used for extraction of total DNA. No specific permissions were required for these locations/ activities. This species was ranked ''Endangered B2ab(iii)'' in ''The IUCN Red List of Threatened Species TM 2013'' based on an outdated information that this species was only found at the type locality (nsc) [31] (http://www.iucnredlist.org). Actually, in addition to the type locality, this species was reported later in 2009 from a variety of localities including the Baiyunshan Mountains (bys), Shirenshan Mountains (srs) and Laojunshan  Mountains (ljs) [32] (Fig. 1), which indicated that this species doesn't meet the criteria of critically endangered or endangered defined in The IUCN Red List of Threatened Species TM 2013, therefore this species should be removed from the red list. However, the new data was not taken into account when this species was ranked ''Endangered B2ab(iii)'' in ''The IUCN Red List of Threatened Species TM 2013''. The specific location (GPS coordinates) of our study was given in Table 1.

Sampling and laboratory protocols
Our sampling covers the entire known distribution of this species. Furthermore, to make an extensive sampling, we explored the whole Tsinling Mountains, and fortunately collected this species at two locations where the distribution of this species has not been recorded. A total of 99 samples were collected from 6 localities during 2011 and 2013 (Table 1, Fig. 1). Eight samples of the alpine toad Scutiger boulengeri were collected from Jone County (34.539922N 103.491647E), southern Gansu Province, China. S. boulengeri will be used as outgroup in phylogenetic relationship analysis.

Determination of generation time
A skeletochronological study of longevity of S. ningshanensis was conducted. The clipped phalanges were stored in 95% ethanol solution. The phalanges were washed in running tap water, then decalcified in 3% nitric acid for 12 to 24 hours. The middiaphyseal region of the phalanges was cross-sectioned at 12-20 mm using a microtome and stained with hematoxylin for two min. Sections were examined under a light microscope and the number of lines of arrested growth (LAGs) was counted [37,38]. The number of LAGs represents the age of the toad.

Phylogenetic analyses
The substitution model selection was implemented in jModelTest 2.1.4 [43], the TIM2+I+G model was selected for all datasets by likelihood ratio tests either under the Akaike Information Criterion (AIC) or under the Bayesian Information Criterion (BIC). Bayesian inference (BI) was used to generate a phylogenetic hypothesis of the DNA haplotypes. BI was performed in MrBayes  generations. Two independent runs each with four simultaneous Monte Carlo Markov chains (MCMC) were conducted. The first 25% of generations were discarded as 'burn-in'. The convergence of chains was confirmed until average standard deviation of split frequency is below 0.01 (0.009889) and the potential scale reduction factor (PSRF) is close to 1.0 for all parameters. In phylogenetic analysis S. boulengeri was used as outgroup.

Population structure analyses
The population comparisons using pairwise difference as distance method, and the partition of genetic diversity within and among populations were analyzed by analysis of molecular variance (AMOVA) [45] using Arlequin3.5.1.2 [42] with 10,000 permutations. Mantel tests [46] were also conducted in Arle-quin3.5.1.2 to assess the significance of isolation by distance (IBD) between populations with 10,000 random permutations on matrices of pairwise population F ST values and the geographical distances. Pair-wise F ST values between populations were estimated using Arlequin3.5.1.2, while geographical distances between populations were calculated online at http://www.gpsvisualizer. com/calculators#distance.
The spatial genetic structure of haplotypes was analyzed using the program SAMOVA1.0 [47] (http://web.unife.it/progetti/ genetica/Isabelle/samova.html) with 1,000 permutations. The number of initial conditions was set to 100 as recommended by Dupanloup et al. [48]. The number K of groups of populations ranged from 2 to 4. The K with the highest Fct represents the best number of groups and the best population configuration.

Historical biogeography
The effective population size (N e ) of each clade (geographical group) for coalescent simulations was converted from Theta using the equation h = N e m with m = 0.65610 28 (per lineage per million years) 66 (generation time of S. ningshanensis). The h-values were estimated using maximum likelihood method in the program Lamarc2.7.5 [48]. Total N e was the sum of the N e for all groups and the proportion of total N e that each group comprised were used to scale the branch width of hypothesized trees (models of population divergence) (Fig. 2) [49][50][51]. Branch widths can be controlled by the Adjust Lineage Widths tool (the horizontal ruler) in the Tree Window in the program Mesquite2.75 [52].
A single-refugium hypothesis of population divergence during the Dali glaciation (the last maximum glaciation in China which occurred about 50 ka before present) was tested using the maximum likelihood estimation implemented in Mesquite2.75 [52] to infer the distributional area of the most recent common ancestor (MRCA) of clades. The single-refugium model hypothesizes that all current geographical populations derived from a single refugium at the end of the Dali glaciation and that all population divergences were concurrent and resulted from the fragmentation of a widely distributed common ancestor's range (Fig. 2). Branch lengths are time in generations based on a six-year generation time in S. ningshanensis. The clades were treated as taxon states and each haplotype as character states. The haplotype trees by coalescence within a diverging population tree (model of population divergence) were simulated, and fit of a haplotype tree to a population tree was calculated to search for population trees that optimize fit of gene trees.

Historical demography
Neutrality tests were calculated in Arlequin3.5.1.2 [41], Fu's F S test [53] and Tajima's D [54] were used to detect evidence of recent demographic expansion within each inferred clade, under which negative value are expected [55]. Inference of population expansion events was performed using a mismatch analysis [56,57] using Arlequin3.5.1.2 with the number of bootstrap replicates set to 10,000 to explore the demographic history of the studied populations. A recent growth is expected to generate a unimodal distribution of pairwise differences between sequences [56]. The validity of the expansion model was tested using the sum of squared deviations (SSD) and Harpending's raggedness index (R) between observed and expected mismatches.
The site-frequency spectrum (for segregating sites) was calculated in DnaSP5. 10.01 [41] to detect the excess of singleton mutations. The null hypothesis of the neutral model (constant population size) was rejected when the allelic frequency spectrum for the entire population revealed an excess of singleton mutations. The excess of singleton mutations could be caused by the expansion [58].
Furthermore, the Bayesian Skyline Plot (BSP) approach was used to estimate the demographic history in BEAST v1.8.0 [59]. The log file was analyzed in Tracer. Strict molecular clock model was selected based on the results from Tracer. A mean mutation rate of 0.65% change per lineage per million years was used [60][61][62][63][64][65][66][67]. Two independent BEAST runs from the same XML file were carried out and then the log output files were combined using LogCombiner. The combined log file was analyzed in Tracer to see if the Trace for each parameter has converged well on a stationary distribution.

Detection of cryptic species
To detect the existence of potential cryptic species, uncorrected p-distances [68,69] for all lineage pairs were calculated in PAUP4.0 [70] from all sequences. The average p-distances of all possible pairs of sequences (every sequence pair contains sequences from different lineages) were calculated.

Generation time
The average longevity of tadpoles was estimated by Lu et al. [33]. The whole longevity equals the longevity of adult plus the longevity of tadpole. The average duration of the tadpole stage of S. ningshanensis is three years [33] and the average duration of the adult stage after metamorphosis is six years. The average life span of S. ningshanensis is nine years. The average generation time (GT) of S. ningshanensis is six years.

Genetic variation
A total of 107 samples were sequenced, including 99 samples of S.ningshanensis and eight of S. boulengeri. A total of 1009 base pairs of cyt b gene and 631 base pairs of COI gene was obtained, 67 haplotypes were identified for the combined gene sequences (cyt b+COI) of S. ningshanensis. Of the combined 1640 nucleotide sites, 181 were variable (polymorphic sites or segregating sites), 121 were parsimony informative, and 60 were singleton variable. The haplotype diversity of total and most sampled population was very high, however, the nucleotide diversity of every sampled population was low (Table 2). Among the total haplotypes, 82.09% are private haplotypes (Table 1).

Phylogenetic relationships among haplotypes
BI analysis revealed three distinct clades (lineages) (hby-lfy, nsc, ljs-bys-srs) in S. ningshanensis. Clades hby-lfy and nsc are well supported with posterior probabilities of 1 and 0.99. On the other hand clade ljs-bys-srs only has a posterior probability of 0.74 (Fig. 3). An important feature of these trees was that the components of each clade showed a strong geographical association. All haplotypes of clade hby-lfy were from western localities, all haplotypes of clade nsc were from locality nsc, all haplotypes of clade ljs-bys-srs were from the eastern localities (Fig. 1).

Population structure and genetic differentiation
AMOVA analysis suggested that majority of the variation occurred among populations (82.15%) while differentiation within populations only contributed 17.85% to the total ( Table 3). The fixation index, a measure of population differentiation due to genetic structure, was highly indicating a significant genetic differentiation among populations (p-value = 0.060.0) ( Table 3).  Population comparisons showed a significant genetic differentiation (F ST ) between most local populations except the population comparison between hby and lfy ( Table 4).
Mantel test results showed a significant correlation between the pairwise calculated genetic distance and pairwise calculated geographical distance (Table 5) of the populations (regression coefficient = 0.001286, correlation coefficient = 0.77051, p (r rand $r obs ) = 0.0185,0.05), indicating the existence of isolationby-distance pattern of genetic divergence for cyt b+COI sequence, which suggests that the distribution of genetic variation is due to geographical separation rather than natural selection. The Mantel test results provided the evidence for large-scale geographical population structure in this species.
Results from SAMOVA analysis indicated that the highest F CT equals to 0.84285 (p (rand.value $ obs.value) = 0.0166260.0) when K = 3, showing that the most likely number of populations is three.

Historical biogeography
The observed value of s is 55 which doesn't fall within the model of single-refugium indicating that the fragmentation model of population divergence was rejected (Fig. 4), that is, the current three lineages were derived from multiple refugia.

Demographic history
The results of mismatch distribution showed that the p-values for SSDs and Rs of the total population and all clades were larger than 0.05, indicating a stable population size in the past; moreover, mismatch distribution in the total population and the clade hby-lfy showed frequency distribution of pairwise difference with multiple peaks (Fig. 5), on the other hand, the p-values for Tajima's D of total population, clade hby-lfy and clade nsc were larger than 0.05 also indicating a stable population size in the past. However, the p-values for Fu's F S s of all clades were smaller than 0.01, and the p-values for Tajima's D of clade ljs-bys-srs was smaller than 0.05, indicating all clades experienced a sudden population expansion (Table 6); on the other hand, mismatch distribution in clades nsc and ljs-bys-srs showed frequency distribution of pairwise difference with single peak (Fig. 5), indicating an sudden population expansion in the past. The two statistics, Tajima's D value and the Fu's F S , are sensitive to   bottleneck effects or population expansion which leads to the more negative values of Tajim's D and Fu's F S [71][72][73][74]. The allelic frequency spectrum for the entire population and all clades revealed an excess of singleton mutations and doesn't fit the neutral model (Fig. 6). The excess of singleton mutations could be caused by the expansion [58].
Bayesian skyline plot (BSP) further estimated the demographic history. The effective sample size (ESS) for all parameters of the BSP was greater than 200, showing that the 20 million generations were sufficient to determine the demographic history for each examined lineage. All lineages and the total population experienced population expansion. The hby-lfy lineage, nsc lineage and the total population experienced quick population growth after about 40,000 years ago, while the ljs-bys-srs lineage experienced a slow population growth after about 15,000 years ago. Noticeably, the hby-lfy lineage and the total population experienced a recent sharp decline in population, the ljs-bys-srs lineage showed a more recent population decline, while the nsc lineage maintained basically constant population size after the population expansion (Fig. 7).

Genetic distances between lingeages
The lineage hby-lfy and the lineage nsc were highly divergent from each other with an uncorrected p-distances of 4.2%, the lineage hby-lfy and ljs-bys-srs were also highly divergent with an uncorrected p-distances of 4.37%. Similarly, there was an high divergence between the lineage nsc and ljs-bys-srs, with an uncorrected p-distances of 2.37%. These values were similar to those p-distances (3%) reported between two different frog species [75]. Therefore S. ningshanensis probably contains three cryptic species or at least three subspecies.

Genetic diversity
Our results showed high levels of genetic diversity in S. ningshanensis. High-level genetic diversity in narrowly endemic species could be associated with the factors such as the long species history, the high-level gene flow, having experienced a recent contraction in population, multiple founder events, the maintenance of genetic diversity in refugia, the hybridization, and the ability to survive in a range of different habitats. [76][77][78]. All these factors except hybridization might contribute to the high level of genetic diversity in S. ningshanensis. This species has probably undergone a long history through which many mutations were accumulated since their origin. High-level gene flow is beneficial to fix the mutations in population. Multiple founder events were another alternative explanation for the high-level genetic diversity in this species. There is a possibility that multiple founder events, which accumulated more and more mutations, occurred in the history of this species. The two congeners of S. ningshanensis are S. boulengeri and S. liupanensis. The nearest localities of these two species are 167 km or 361 km respectively away from the locality hby of S. ningshanensis, which makes it impossible for the occurring of hybridization among these three species [79,80]. Therefore, hybridization may not contribute to the high-level genetic diversity in S. ningshanensis. The distributing region of S. ningshanensis include a range of different habitats such as high mountains with different elevation, streams, and different vegetations. High-level diversity in habitats also contributed to the high-level genetic diversity.

Population structure
Significant population structure occured based on the statistics pairwise differentiation. Most pairwise F ST values are high and Table 6. Mismatch distribution analyses and neutrality test of S. ningshanensis. statistically significant. Thus, the Ningshan alpine toad appears to exhibit substantial population differentiation across the Tsinling Mountains.
Mantel test results showed a significant correlation between the genetic distance and geographical distance of the populations, indicating the presence of IBD pattern of genetic divergence for cyt b+COI sequences, suggesting that the distribution of genetic variation is due to geographical separation rather than natural selection. The Mantel test results provided the clear evidence for large-scale geographical population structure in S. ningshanensis. It is not possible that a significant Mantel test provided the evidence for discontinuity in the distribution of genetic variation. It rather showed a continuous distribution of the variation due to individuals mating preferentially with individuals from nearby populations [81].
AMOVA results indicated that 82.15% genetic variation occurred among populations, while differentiation within the populations only made 17.85% contributions. The high genetic variation among populations affirmed the presence of phylogeographic structure in S. ningshanensis. SAMOVA results and the phylogenetic relationship analysis further affirmed the existence of phylogeographic structure in this species.
Topography and Pleistocene climate changes drive population genetic differentiation forming genetic structure pattern [90][91][92]. East Asia including China has undergone a series of cooler-drier climate cycles in the last 15 million years [93]. Dramatic climatic changes have caused the extinction of many organisms [94]. S. ningshanensis that distributed in the areas with low elevation disappeared during the interglacial in the Quaternary since it is an alpine species. In addition, S. ningshanensis retreated to a few refugia during glacial period. At the end of the Dali glaciation, S. ningshanensis experienced a rapid population expansion which occurred about 40,000 years before present, and it is experiencing a population contraction now probably due to anthropogenic disturbance. Topographic features, climatic fluctuation and anthropogenic activity together led to the current patchy geographic pattern for S. ningshanensis. This geographical distribution pattern also seen in other montane organisms [95,96].
Finally, there is the possibility that unsampled populations between nsc and ljs may be genetically intermediate among the three groups.

Historical biogeography
The uplifts of the Tsinling Mountains promoted genetic differentiation among lineages of S. ningshanensis. The Tibetan Plateau experienced three phases of rapid uplifts [97,98]. The uplifts of the Tibetan Plateau impacted the environment of the surrounding areas including the Tsinling Mountains [97]. The Tsinling Mountains experienced similar uplift process, and three phases of uplifts occurred in Pliocene, Early Pleistocene and Holocene respectively [99][100][101][102]. Coalescent simulations indicated multiple refugia in S. ningshanensis though we could not determine the number and location of the refugia based on the statistical phylogeography analyses. Most species of the genus Scutiger inhabit in the eastern Tibetan Plateau [81], we guess that S. ningshanensis evolved when the Scutiger species dispersed eastward along the Tsinling Mountains. The hby-lfy lineage split apart at first, split between the lineages nsc and ljs-bys-srs occurred when this species continued dispersing from west to east on the Tsinling Mountains, while no further split occurred after the second split.

Demographic history
Most demographic analyses revealed a sharp population expansion in all lineages of S. ningshanensis, and the expansion began simultaneously 40,000 years before present, which corresponds to the end of the Dali glaciation. Retreat of glacier led to population expansion in species on the Tsinling Mountains. A noteworthy phenomena is that S. ningshanensis is experiencing a distinct recent population contraction (though the lineage nsc is exceptional) as shown by the BSP analyses (Fig. 7), there is a possibility that anthropogenic disturbance resulted in the contraction in population size of this species. Multiple uplifts of the Tsinling Mountains and fluctuations in population size of S. ningshanensis associated with glacial-interglacial cycles led to increases or decreases in the levels of genetic variation and coalescence times [23,25,28].
Putative cryptic species S. ningshanensis probably contains at least two cryptic species or subspecies based on the p-distance analysis though there is not difference in morphology among these cryptic species or subspecies. The geographical distances between the local populations are long enough for occurrence of gene flow break in S. ningshanensis since it is philopatric to breeding sites. In addition, high peaks and deep valleys also contributed to the break of gene flow between populations. Subsequently, poor level of gene flow led to speciation events. However, there are limitations in taxonomic consequences based on only one taxonomic discipline. We will further confirm the cryptic speciation using data from nuclear DNA and ecological niche modeling (ENM) in the future.

Conservation and management implications
All potential cryptic species should be considered for conservation. Different conservation strategies should be accepted for different species, therefore, it is inappropriate to protect a cryptic species complex using a single conservation strategy [103]. It is indispensable to understand and quantify biodiversity so that we can better explain and at last carry out conservation [103]. The distribution area of this species is limited to a narrow zone along the Tsinling Mountains, if this species is a cryptic species complex itself, the distribution area of each cryptic species is much smaller which increases the risk for extinction. Effective protection measures are to be provided and carried out. Conservation should be considered for every cryptic species.
Management units (MUs) are used to make conservation strategy for threatened species [104,105]. The three lineages (hby-lfy, nsc, ljs-bys-srs) derived from our study could be used as three different MUs, each MU has its own unique genetic structure. Every MU should be considered when conservation policy is made.

Conclusions
Three lineages were detected. The phylogenetic relationship analyses showed significant phylogeographic structure. The population expansion or contraction and genetic differentiation between populations or lineages could be explained by topography and the repetitive uplifts of the Tsinling Mountains and the climatic cycles during the late Pliocene and Pleistocene. S. ningshanensis experienced a rapid population expansion about 40,000 years before present. The current decline in population size was probably caused by anthropogenic disturbance. Current populations of S. ningshanensis are from different refugia though the location of these refugia could not be determined in our study. Topography is very important in shaping genetic connectivity. Topography, climatic changes and repetitive population expansion/contraction together led to the high level of genetic variation in S. ningshanensis. At least two putative cryptic species were detected. A total of three MUs were determined, which must be considered when conservation policy is made in future.