Molecular Phylogeography and Population Genetic Structure of an Endangered Species Pachyhynobius shangchengensis (hynobiid Salamander) in a Fragmented Habitat of Southeastern China

The salamander Pachyhynobius shangchengensis (Hynobiidae) is a vulnerable species restricted to a patchy distribution associated with small mountain streams surrounded by forested slopes in the Mount Dabieshan region in southeastern China. However, molecular phylogeography and population genetic structure of P. shangchengensis remain poorly investigated. In this study, we explored the genetic structure and phylogeography of P. shangchengensis based on partial sequences of the mitochondrial DNA (mtDNA) cytochrome b and cytochrome c oxidase subunit I genes. Fifty-one haplotypes and four clades were found among 93 samples. Phylogenetic analyses revealed four deeply divergent and reciprocally monophyletic mtDNA lineages that approximately correspond to four geographic regions separated by complicated topography and long distances. The distinct geographic distributions of all lineages and the estimated divergence time suggest spatial and temporal separation coinciding with climatic changes during the Pleistocene. Analysis of molecular variance indicated that most of the observed genetic variation occurred among the four groups, implying long-term interruption of gene flow, and the possible separation of P. shangchengensis into four management units for conservation.


Introduction
Molecular ecology primarily aims to understand the influence of abiotic factors, such as altitude, topography, and glacial history, on the spatial distribution of genetic variations [1]. Recently, researchers have studied the effect of landscape variables, such as topography and altitude, on the geographical distribution of genetic variation in the emerging field of landscape genetics [2]. Landscape characteristics can affect the proportion of suitable habitats, migration patterns, and the genetic divergence of populations [2,3]. For terrestrial species, especially for amphibian, species diversification drived by landscape variables can occur by the formation of sky islands in which hot, dry, deep valleys serve as barriers to gene flow [4], as well as the height of mountains forming a barrier to dispersal for amphibian that live in the valleys [5]. Mountain ridges have also been shown to be an important barrier to amphibian dispersal and gene flow [2,6]. Therefore, a complex, microhabitat-rich topography, could effect genetic diversity and phylogeographic structure of animal habitats in these areas [7].
Using phylogeographical tools to analyze the effects of landscape characteristics on species distributions over large spatial scales has provided remarkable insight into the spatial patterns of genetic diversity [8]. Using highly variable genetic markers and a dense sampling regime across a small, topographically diverse region enables investigation of the localized effects of geography on genetic diversity and connectivity across the landscape [1].
Besides landscape characteristic, climatic changes have also caused montane species to expand, change, or be in contact with each other along latitudinal or elevational gradients associated with Pleistocene glacial cycles [9][10][11]. "East Asia is characterized by a mosaic of mountains and likely experienced a relatively mid-Pleistocene climate" [12]. The Dabieshan Mountains are connected to the eastern end of the Tsinling Mountains by the Tongbaishan Mountains, and are located in the eastern part of China. Though most mountains in this area were not glaciated during the Pleistocene [13,14], this region experienced climatic fluctuations which probably impacted species distributions, demography and diversification [15]. Phylogeography and population genetic structure of P. shangchengensis which lives in here should be affected by climatic changes.
Shangcheng Stout Salamander P. shangchengensis (Hynobiidae) is an endemic species in china, its distributions restrict in the Mount Dabieshan regions, which was described by Fei et al. [16] from Mount Huangbaishan, Shangcheng County, Henan Province (holotype) and Mount Jingangtai, Shangcheng County (paratypes). It can be found in patchy habitat on the Dabieshan Mountains, southeastern China [17,18]. Pachyhynobius shangchengensis has low vagility, its habitat is separated by valleys and low lands, and its distributions are getting smaller and smaller [17]. Chen et al. [19] described Hynobius yunanicus based on specimens from Huangbaishan, Shangcheng County. Hynobius yunanicus differs from P. shangchengensis mainly in having little white spots on deep brown dorsal side, in lacking premaxillary fontanelle on the skull, and in lacking connection between maxillary and pterygoid. However, evidences from karyotypic and phylogenetic analysis rejected the validity of H. yunanicus [20]. Therefore, Hynobius yunanicus is a synonym of P. shangchengensis, and our sampling sites should include Huangbaishan.
In this study, we test whether the topography of the mountains affected the population genetic structure of the P. shangchengensis according to partial sequences of the mitochondrial DNA cytochrome b (mtDNA cyt b) and cytochrome c oxidase subunit I (mtDNA COI) genes. We also determined whether climatic oscillations during glacial periods in the Quaternary affected the distribution of P. shangchengensis.

Sampling, DNA extraction, polymerase chain reaction, and sequencing
This study was approved by the Institutional Animal Care and Use Committee (IACUC) of Shaanxi Normal University and Chinese Academy of Sciences. A total of 93 P. shangchengensis individuals were collected from seven locations on the Dabieshan Mountains, Southeast China (Table 1; Figure 1). Samples were obtained upon capture by toe  Samples were also permitted by following authorities: management committee of Jingangtai Nature Protection Area; management committee of Shangcheng Stout Salamander Nature Protection Area (Henan Shangcheng); management committee of Tiantangzhai National Forest Park; management committee of Yaoluoping Nature Protection Area.
Tissue samples were preserved in 95% ethanol and stored at -20 °C. Total genomic DNA was extracted through a standard phenol: chloroform method [21]. A continuous fragment (942 bp) of the mitochondrial cytochrome b gene was amplified by polymerase chain reaction (PCR) (MyCycler Thermal Cycler), with outer primers L 1 , H 1 and inner primers (nested primer) L 2 , H 2 ( Table 2). A continuous fragment (1011bp) of the mitochondrial cytochrome c oxidase subunit I (COI) was amplified by PCR (MyCycler Thermal Cycler), with forward primer FZ and reverse primer RZ ( Table 2). The PCR products were checked in a 1% agarose gel and purified through a TIANquick Midi Purification Kit (Tiangen, Beijing, China) according to protocol recommendations. Sequencing reactions were performed with the PCR primers through an ABI Prism BigDye Terminator Cycle Sequencing-Ready Reaction Kit on an ABI 3730XL sequencer. All sequences were deposited in the GenBank database under accession numbers KC162002-KC162083 (COI) and KC162084-KC162162 (cyt b).

Nucleotide polymorphism
The sequences were aligned with Clustal X1.83 [22]. The aligned sequences were edited using the program BioEdit 7.0.9.0 [23]. Haplotype inference was conducted through Collapse 1.2 (http://darwin.uvigo.es). To examine whether the two analyzed regions (cyt b and COI) can be combined into a larger data matrix [24], we performed a partition-homogeneity test using 1000 replicates as implemented in PAUP 4.10b [25]. The combined data were further analyzed because the result of the partition-homogeneity test was not significant.
The number of variable and parsimony-informative sites was determined using the program DnaSP 5.10.01 [26], and haplotype diversity (H d ) and nucleotide diversity (P i ) were determined through Arlequin 3.5.1.2 [27].

Phylogenetic structure
The phylogenetic relationship among haplotypes was estimated through maximum likelihood (ML) analyses in PAUP*4.0b10 [25], as well as Bayesian analyses in MrBayes 3.0 [28] with 3,000,000 generations. For the Bayesian analyses, MrModelTest 2 [29] was used to find the best-fit substitution model, and GTR+I+G model was performed. For the maximum likelihood analyses, MODELTEST [30] was used to find the best-fit substitution model, and GTR+I+G model was performed. The confidence level of ML trees was accessed by 1000 bootstrap replications. Hynobius chinensis and Hynobius guabangshanensis were used as the outgroup. The COI and cyt b sequences of H. chinensis were downloaded from GenBank with accession number HM036353, and the COI and cyt b sequences of H. guabangshanensis were downloaded from GenBank with accession number EF616473 (cyt b) and FJ913877 (COI). We also used NETWORK 4.5.0.2 [31] to draw a medianjoining network to analyze the relationships among detected haplotypes.

Analyses of geographic structuring
The population and phylogroup comparisons using pairwise difference and the partitions of genetic diversity within and among populations were analyzed through analysis of molecular variance (AMOVA) [32] using Arlequin 3.5.1.2 [27] with 10,000 permutations.
The spatial genetic structure of haplotypes was analyzed through SAMOVA 1.0 [33] (http://web.unife.it/progetti/genetica/ Isabelle/samova.html) with 1000 permutations. The number of initial conditions was set to 100 as recommended by Dupanloup et al. [33]. The number K of groups of populations was set to from 2 to 6 respectively. The K with the highest FCT represents the best number of groups and the best population conFigureuration. This program implements an approach to define groups of populations that are geographically homogeneous and maximally differentiated from each other. The method is based on a simulated annealing procedure that aims to maximize the proportion of total genetic variance caused by differences between groups of populations (FCT).

Divergence time estimate
The approximate divergence times were estimated for the lineages for P. shangchengensis in BEAST 1.6.1 [34]. Except those for outgroups, all haplotype sequences were used in the analysis. A Bayesian Markov chain Monte Carlo approach with an uncorrelated log-normal relax molecular clock was used in BEAST 1.6.1 [34]. Two independent runs were performed, each of which was composed of 120 million generations, with sampling every 1000 generations. A burn-in was set to 10% of the samples. To check for stationarity, the results were displayed in TRACER version 1.5 [35]. LogCombiner 1.4.7 [36] was used to combine both runs. TreeAnnotator 1.4.7 [36] was used to annotate tree information, and FigureTree 1.1.1 [37] to visualize tree information.
A molecular evolutionary rate of the mitochondrial genome for hynobiids (0.64% per Myr per lineage) was proposed by Weisrock et al. [38]. This evolutionary rate was frequently used for hynobiids mitochondrial DNA data [39][40][41][42]. Thus, we used 0.64% per Myr per lineage to estimate the divergence between any major clades.

Pattern of isolation by distance
Mantel tests [43] were conducted in Arlequin 3.5.1.2 [27] to assess the significance of isolation by distance between populations with 5000 random permutations on matrices of pairwise population F ST and the geographical distances. Pairwise F ST values between populations were estimated through Arlequin 3.5.1.2 [27], whereas straight line geographical distances between populations were calculated online at http://www.gpsvisualizer.com/calculators.

Demographic history analysis
We applied Neutrality tests through the program Arlequin 3.5.1.2 [27] as an assessment of possible population expansion. Under the assumption of neutrality, a population expansion produces a large negative value of Fu's F S test [44] and Tajima's D [45]. Tajima's D and Fu's F S are sensitive to bottleneck effects or population expansion, causing these values to be more significantly negative [46][47][48][49]. Fu's F S is particularly sensitive to recent population growth [44].
Population expansion events were determined through mismatch analysis [50] using Arlequin 3.5.1.2 [27], with the number of bootstrap replicates set to 5000 to explore the demographic history of the studied populations. The parameters of demographic expansion were also estimated. Recent growth is expected to generate a unimodal distribution of pairwise differences between sequences [50]. The validity of the expansion model was tested by using the sum of squared deviations (SSD) and Harpending's raggedness index (R) between observed and expected mismatches. The formula T=τ/ (2ut) was used to estimate the time of the population expansions [50] based on the generation time (3.5 years) [51], t is the date of growth or decline (mutational time), τ is the mode of mismatch distribution (evolutionary time), and u is the mutation rate per sequence and per generation [52].
The Bayesian skyline plot (BSP) was used to estimate the demographic history of P. shangchengensis using the program BEAST 1.6.1 [34]. A piecewise-constant skyline model was selected, and a relaxed uncorrelated log-normal molecular clock was used with the mutation rate of 0.64%/MY for P. shangchengensis as suggested by Weisrock et al. [38]. Tracer 1.5 was used to reconstruct the demographic history through time.

Genetic variation
The total number of sites (excluding sites with gaps/missing data) was 1953, of which 1011 bp were sequenced for the COI gene and 942 bp for the cyt b gene. A total of 195 polymorphic sites were found, of which 182 were parsimony-informative and 13 were singleton-variable. These polymorphic sites identified 51 haplotypes within 93 individuals from seven localities (Table 1; Figure 1). Each sampled population and the total population have high H d , accompanied by very low P i (Table 1).
For the spatial AMOVA, with K increased from 2 to 6, the F CT value was highest (F CT = 0.9421) when K = 4. Thus, the SAMOVA tests revealed the number of significant phylogeographic groups (K = 4) .

Pattern of isolation by distance
Mantel test results showed significant correlation between the pairwise calculated genetic distance and pairwise calculated straight line geographical distance of the populations (correlation coefficient = 0.6406, p < 0.001), indicating the presence of isolation-by-distance. This finding suggests that the distribution of genetic variation is due to geographical separation. The Mantel test results provided evidence for largescale geographical population structure in this species.

Population Comparisons
Population comparisons showed more or less significant genetic differentiation (FST) between most local populations (Table 4). Thus, a long-term interruption of gene flow among all clades was also evidenced by the relatively high FST values.

Demographic inferences and divergence time
The results of neutral test analyses of clade A, B, C indicated that both Tajima's D and Fu's Fs were negative, and some values have highly significant, except clade D possessed a positive value (Table 5). Mismatch distribution analyses showed a unimodal frequency distribution of pairwise differences in clade A, clade B and clade C ( Figure 5). All above results suggest demographic expansion. The estimated expansion time of above clades was 0.06-0.03 Myr in the Late Pleistocene, the results were consistent with the analysis of the BSP (Table 5; Figure 6). And a sudden expansion was identified between 0.05-0.0025 Myr by BSP ( Figure 6). However, both mismatch distribution analyses and the neutrality tests rejected a sudden population expansion in the clade D and total population (Table 5; Figure 5).
The results of analyses in the program BEAST inferred that the estimated age of the origin of P. shangchengensis on Mount Dabieshan in China to be 23.39 Myr. The divergence time between clade A and clade B was calculated to have taken place in the early Pleistocene 2.11 Myr, with a 95% highest posterior density (HPD) of 1.54-2.74 Myr. The divergence time between clade C and clade A+B was

Pre-Pleistocene split and geologic history
The various drivers of species divergence associated with topography seem to play roles in the evolution of P. shangchengensis. In our study, we can see that phylogenetic analyses support four major clades. The median-joining network yielded four unconnected subnetworks corresponding to the four clades in the phylogenetic tree, and there are no shared haplotypes between clades (Figures 2, 3, 4 ). Such distributions of mitochondrial DNA haplotypes of P. shangchengensis may be interpreted as being the result of population isolation because of their specific biological habits.
There are several reasons for Amphibians which are particularly sensitive to effects of topographic and altitudinal variation. Amphibians are generally highly site philopatric and poor dispersers [53][54][55][56]. Because of dessication and predation risks associated with terrestrial dispersal [57][58][59] and slow terrestrial locomotion [59], low vagility in amphibians is often attributed to dependence on moist habitats or wetland corridors for dispersal [59,60]. Thus, because of complex topography in the Mount Dabieshan regions, consequent range expansion and dispersal of individuals away from their natal sites are generally expected to be limited in P. shangchengensis. And this is the reason that P. shangchengensis can be found in patchy habitat on the Dabieshan Mountains [17,18,61].
Topography can drive divergence patterns [62], Pleistocene climatic fluctuations associated with cyclical glaciation events can drive divergence patterns too, even in lower latitudes [63,64]. China and its neighboring areas in East Asia have experienced a development of cooler and drier climates within the last 15 Myr, although most of China has never been covered by ice sheets [65]. Furthermore, tremendous climatic changes, particularly the Quaternary glaciations, have made many plants and animals extinct and influenced the evolution and distribution of many plants and animals in China and its neighboring areas, particularly during the Quaternary [66]. In our study, BEAST analysis showed that the estimated age of P. shangchengensis was 23.39 Myr (in the early Miocene), between the main clades, the divergence time taken place in the Pliocene, and by Bayesian skyline plot (BSP) a sudden expansion occurred in the Late Pleistocene. Therefore, climatic fluctuations probably impacted the distribution, demography and diversification of the species. We infer that P. shangchengensis which lived in the low lands disappeared during the interglacial in the Quaternary, because it could not adapt to hot and dry climate. By contrast, Pachyhynobius shangchengensis, which lived in higher elevations, survived because these areas were suitable for survival. Thus, the populations of P. shangchengensis have been isolated in the fragmented mountain habitats in the last interglacial to the present. This pattern also existed in other montane organisms [67,68].

Conservation and management implications
The two goals of any conservation program are to maintain the genetic diversity of species for long-term evolutionary success and ensure their survival [20,69]. The number of P.
Management units (MUs) are commonly used designations for threatened or endangered taxa [70,71]. The data of our study can be used to establish MUs because these units are defined by either reciprocal monophyly in mtDNA or substantial allele frequency divergence at nuclear loci [70]. Thus, the four populations of P. shangchengensis can be considered MUs because genotypes in the four populations are closely related but not shared. Any conservation policy should concentrate on protecting the distinct populations with similar MUs similar to the conservation efforts for the Tibetan gazelle in China [15]. We would do our best to recommend to protect P. shangchengensis and avoid its extinction, we should take necessary measures, such as forbidding to kill any individual species in distribution areas of P. shangchengensis.

Conclusions
In conclusion, Pachyhynobius shangchengensis has significant phylogeographic structure, the topography of the Dabieshan Mountains significantly affects the population genetic structure of it, and climatic oscillations during glacial periods in the Quaternary affected the distribution of this species.