Skip to main content
Download PDF

Long-term sky islands generate highly divergent lineages of a narrowly distributed stream salamander (Pachyhynobius shangchengensis) in mid-latitude mountains of East Asia

BMC Evolutionary Biology volume 19, Article number: 1 (2019) Cite this article

Abstract

Background

Climate oscillation may have a profound effect on species distributions, gene flow patterns and population demography. In response to environmental change, those species restricted to montane habitats experienced expansions and contractions along elevation gradients, which can drive differentiation among sky islands.

Results

The Shangcheng stout salamander (Pachyhynobius shangchengensis) is a cool stream amphibian restricted to high-elevation areas in the Dabie Mountains, East China. In the present study, we used mtDNA genes (Cyt b and ND2) of 193 individuals and 12 nuclear microsatellite loci genotyped on 370 individuals, representing 6 populations (JTX, KHJ, MW, TTZ, BYM and KJY) across the taxon’s distribution area, to investigate their genetic variation and evolutionary history of P. shangchengensis. Most populations showed unusually high levels of genetic diversity. Phylogenetic analyses revealed five monophyletic clades with divergence times ranging from 3.96 to 1.4 Mya. Accordingly, significant genetic differentiation was present between these populations. Bayesian skyline plot analyses provided that all populations underwent long-term population expansions since the last inter-glacial (0.13 Mya ~ 0.12 Mya). Msvar analyses found recent signals of population decline for two northern populations (JTX and KHJ) reflecting a strong bottleneck (approximately 15-fold decrease) during the mid-Holocene (about 6000 years ago). Ecological niche modelling has shown a discontinuity in suitable habitats for P. shangchengensis under different historical climatic conditions.

Conclusions

Our results suggest that the niche conservatism of P. shangchengensis and sky island effects may have led to long-term isolation between populations. In sky island refuges, the mid-latitude Dabie Mountains have provided a long-term stable environment for P. shangchengensis, which has led to the accumulation of genetic diversity and has promoted genetic divergence.

Background

Geological, climatic, and hydrological histories are the major drivers that have determined present-day patterns of biodiversity [1]. These factors can change the diversification at different temporal and geographical scales [2, 3]. For example, the uplifting of the Qinghai-Tibet Plateau caused habitat fragmentation, limitations to dispersal and the formation of gene flow barriers that accelerated divergence and even led to speciation [4, 5]. During the last million years, climatic oscillations, such as the glacial expansions and retreats, may have led to current patterns of species distributions and inter-population gene flow [6,7,8]. Hydrological histories might change dramatically under the influence of climatic oscillations and geological changes, which can influence the phylogeographic patterns drastically by the barrier effect, especially for aquatic organisms [9, 10]. In addition, during the Quaternary period, severe climatic oscillations made a profound impact on the phylogeography of species at mid-latitudes, especially for those endemic to mountainous areas [11, 12]. Therefore, the effects of geological, climatic, and hydrological histories were perhaps most pronounced in mid-latitude mountains where species adapted to local environmental conditions may have experienced alternating periods of isolation and connection during climatic fluctuations [10, 13,14,15]. The genetic consequences of these historical events are manifested in the phylogeographic patterns and levels of genetic variation across multiple taxa [6,7,8].

Sky islands are montane regions isolated from one another by intervening valleys with drastically different environmental conditions, in a similar way to how the sea separates oceanic islands [16]. In sky islands, changing climatic conditions can cause favorable environments for species or populations to expand or contract along an elevational gradient, leading to alternating periods of isolation and connection [6, 10,11,12,13,14, 17]. Species restricted to sky islands, especially aquatic taxa, often show unique patterns of population genetic structure [10, 13, 14, 17] influenced by historical climate-induced distributional shifts, and are highly vulnerable to future climate changes [6, 18]. Therefore, montane habitats often become hotspots of endemism, intraspecific genetic diversity and species richness, and provide the opportunity to study how different evolutionary processes lead to diversification [10, 12,13,14, 18, 19].

The Dabie Mountains, located at mid-latitude in East Asia, are composed of a chain of ancient, isolated, low-middle elevation massifs (Fig. 1). The Dabie Mountains are located in the ecotone of subtropical evergreen broad-leaved forest and the warm-temperate deciduous broad-leaved forest zone. The climate in this region is transitional from subtropical to temperate, with annual average temperature of 12.5 °C and mean annual precipitation of 1832.8 mm [20]. The Shangcheng stout salamander (Pachyhynobius shangchengensis) is a stream salamander narrowly distributed in the Dabie Mountains and is endemic to the cool and oxygen-rich streams above 500 m [21]. Presently, P. shangchengensis is found distributed as six isolated areas/populations (Jiaoyuan-Tanghui-Xiaolongtan, JTX; Kangwangzhai- Huangbaishan- Jiufengjian, KHJ; Mazongling-Wochuan, MW; Tiantangzhai, TTZ; Baimajian-Yaoluoping-Mingtangshan, BYM; Kujingyuan, KJY) (Fig.1). In previous studies, distinct genetic divergence was found between some of these areas, however, these studies differed with respect to the mitochondrial gene and areas analysed (i.e. control region and JTX, KHJ, MW, TTZ, BYM [22], and Cyt b, COI and JTX, KHJ, TTZ, BYM [23]). In the mid-latitude areas, especially in mountains, species adapted to the local environment often show unique phylogeographic patterns, such as North American Plethodon [10, 13, 14]. Based on the above considerations, we have argued that a similar evolutionary course may also exist in the mid-latitude mountain areas of East Asia. In view of the distribution area characteristics, together with existing phylogeographic results for P. shangchengensis, we predicted that P. shangchengensis may be an ideal model organism to understand how climate change can impact genetic variation and phylogeographic patterns of species in mid-latitude mountains of East Asia.

Fig. 1
figure 1

Sampling area and populations of P. shangchengensis. Sample sites are represented as triangles or dots in different colors. The map derived from https://glovis.usgs.gov/. The basic vector layers of rivers derived from http://www.resdc.cn/

Full size image

Here we sampled P. shangchengensis throughout its range and used multiple molecular markers combined with GIS–based environmental niche analyses to better understand population genetic relationships and evaluate diversification processes. Firstly, we used statistical phylogenetic methods to test whether each geographically isolated sky island population is reciprocally monophyletic and whether the patterns of phylogenetic diversification identified are consistent with past climatic oscillation. Secondly, we examined whether the isolated sky island populations maintained genetic diversity effectively. Lastly, we used ecological niche modelling (ENM), combined with paleoclimate data to infer historical causes for the isolated distribution pattern of P. shangchengensis.

Methods

Sampling collection

Samples were collected from 370 individuals from 6 populations (JTX, KHJ, MW, TTZ, BYM and KJY), throughout the distribution range of P. shangchengensis between 2012 and 2015 (Fig. 1). Sampling involved capturing adults with dip nets and removing the tip of the tail (about 1 cm) prior to releasing the individuals. Samples were preserved in absolute ethanol in the field, and then stored at − 80 °C until used. For each population, we sampled no less than 20 individuals (Table 1). Individuals caught twice in different years (identified by their docked tails) were not sampled.

Table 1 Descriptive statistics for mitochondrial genes and 12 micorsatellite loci for P. shangchengensis populations in Dabie Mountains
Full size table

DNA extraction, sequencing, and microsatellite genotyping

Total DNA was extracted from samples using a standard proteinase K/phenol-chloroform protocol [22], and purified by EasyPure PCR Purification Kit (TransGene). Cytochrome b (Cyt b) and NADH dehydrogenase 2 (ND2) were amplified using primers designed from Ranodon sibiricus (NC004021) and P. shangchengensis (NC008080) mitogenome using Primer v5 [23] (Additional file 1: Table S1). PCR reaction mixtures (50 μL) consisted of 2 μL genomic DNA (concentration 10–50 ng/μL), 5 μL 10× buffer, 2 μL of 2.5 mM MgSO4, 4 μL of 2 mM dNTPs, 2 U Taq polymerase (TransGen Biotech, Beijing), 0.6 mM of each primer and sufficient pure molecular biology grade water. The amplification program included an initial denaturation step of 95 °C for 5 min followed by 32 cycles of denaturation at 95 °C for 30 s, primer annealing at 55 °C for 30 s and an extension at 72 °C for 80 s, with a final extension at 72 °C for 10 min. PCR products were purified using EasyPure PCR Purification Kit (TransGene), and sequenced on an ABI Prism 3730 automated sequencer using the BigDye Terminator v3.0 Ready Reaction Cycle Sequencing Kit (Applied Biosystems).

Twelve microsatellite loci were selected from previous studies (in Additional file 1: Table S1) [24]. PCR reaction mixtures (25 μL) consisted of 1 μL genomic DNA (concentration 10–50 ng/μL), 2 μL 10× buffer, 1 μL of 2.5 mM MgSO4, 2 μL of 2 mM dNTPs, 1 U Taq polymerase, 0.3 mM of each primer (forward primer fluorescently labeled with FAM, HEX or TAMRA) and sufficient water. The amplification program was conducted with following conditions: 5 min at 95 °C; followed by 35 cycles of 30 s at 95 °C, 20 s at the annealing temperature (Additional file 1: Table S1), 30 s at 72 °C; and 5 min at 72 °C. PCR products were visualized on an ABI 3730 semiautomated sequencer (PE Applied Biosystems) with the GS500 marker as size standard and then analyzed using GeneMarker 1.85 (version 1.3, SoftGenetics LLC) [25].

Data analysis

Sequences were assembled using SeqMan in DNAStar [26]. All DNA sequences were aligned using Clustal X v.2.0 [27], trimmed, and deposited into GenBank (MF581931–MF582316). We estimated population nucleotide diversity (π) and haplotype diversity (Hd) in DnaSP v.5.10. Micro-Checker [28] was used to check the presence of null alleles and genotyping errors in microsatellite data. Deviation from Hardy-Weinberg equilibrium (HWE) for each locus and for all loci for each population was evaluated using exact probability tests implemented in GenePop v.4.2.1 [29]. Several population genetic summary statistics to describe genetic variation were estimated by GenAlEx v.6.5 [30] and GENETIX v.4.02 [31], including mean number of alleles per locus (Na), observed heterozygosities (HO), expected heterozygosities (HE), inbreeding coefficients (F) and pairwise FST.

Phylogeography and population structure

We reconstructed the phylogenetic relationships and divergence time of combined sequences (Cyt b and ND2) using BEAST v.1.8 [32], with R. sibiricus (NC004021) as outgroup. The best partition strategy and nucleotide substitution model for each data partition were determined with Partitionfinder v.2.1.1 [33]. The divergence times were estimated using a Bayesian Markov Chain Monte Carlo (MCMC) method with a strict clock method as implemented in BEAST. Due to the lack of suitable fossil calibrations for our analyses, we placed a lognormally distributed prior (centered at 40.2 Mya and with a 95% HPD between 34.5 and 46.2 Mya) on the age of the node marking the divergence between P. shangchengensis and the outgroup, based on the time calibrated molecular phylogenetic tree of Hynobiidae [34]. We simultaneously estimated topology and divergence dates by performing two independent runs of 20 million MCMC iterations sampling every 2000th iteration with the initial 25% of the samples discarded as burn-in. An uncorrelated lognormal model of lineage variation with an expansion growth tree prior and a strict molecular clock was used [32] based on the model comparisons carried out in Tracer v.1.6 [35] including a comparison with the alternative relaxed molecular clock model. Convergence of the two independent MCMC runs was assessed in Tracer, as were convergence of model parameter values (effective sample size [ESS]) to ensure ESS values > 200. The tree and posterior distribution were summarized with TreeAnnotator v.1.8 [36] and visualized in FigTree v.1.4.3 [37].

A minimum-spanning network based on mitochondrial genes was created with Network v.4.6.1 (Fluxus Technology, Suffolk, UK), using the full median-joining algorithm [38] to visualize relationships of the mitochondrial DNA (mtDNA) haplotypes among localities. Genetic differentiation (FST) between populations was estimated in Arlequin v.3.5.1.3 [39] based on haplotype frequency differences and with 10,000 permutations to assess significance. Homogeneity of haplotype distributions among localities was tested using a single-level AMOVA in Arlequin.

A Bayesian analysis of population structure using STRUCTURE [40] was carried out to estimate the number of potential clusters present in the microsatellite data and to assign individuals to inferred clusters. Ten independent runs were carried for different values of K (the number of clusters) between 1 and 8, using no prior information about individual location, and assuming admixture and correlated allele frequencies. The MCMC analyses were run for a total of 1 million generations discarding the first 100,000 as burn-in. The most likely K explaining the variation in the data was selected by choosing the K result for which the log likelihood [Ln Pr(X/K)] of the posterior probability of the data was the highest [40], and the ΔK statistic [41]. The results were graphically displayed with the software DISTRUCT [42].

A spatial cluster analysis was conducted with the Geneland package [43] in R 3.3.1 (R Development Core Team 2016). The MCMC iterations were set to 100,000 steps sampling every 100th iteration and discarding the first 10% iterations as burn-in period. K values between 1 to 6 were tested with individuals assigned to one K cluster (1 ≤ K ≤ 6) based on their multilocus genotypes and spatial coordinates. To confirm that Geneland was long enough, we carried out 10 different runs and compared the parameter estimates between them (K, individual population membership, maps). The best result was chosen based on the highest average log posterior density.

Population demography

Mismatch distributions were estimated to test the null hypothesis of recent population growth using Arlequin. A goodness-of-fit test was used to determine the smoothness of the observed mismatch distribution (using Harpending’s raggedness index, Rag) and the degree of fit between the observed and simulated data (using the sum of squares deviation, SSD) [44]. Due to their sensitivity to demographic changes, we also estimated Tajima’s D [45] and Fu’s Fs [46] using 10,000 coalescent simulations to assess significance in Arlequin. We also examined past population dynamics of all phylogeographic lineages of P. shangchengensis using the Bayesian skyline plots (BSP) with the model selected by Partitionfinder. BSP were performed in BEAST to reconstruct historical population size dynamics by using the same parameter sets as previous divergence time estimation. Each dataset was run four times for 20 million MCMC iterations, sampling every 2000 iterations and discarding the first 25% of the iterations as burn-in. Tracer was used to check the convergence of the MCMC analyses (ESS values > 200).

The software Msvar v.1.3 [47] was used to characterize the recent demographic history of the six populations of P. shangchengensis using the microsatellite data. The method assumes that a current population (of size N0) passed through a demographic change (a bottleneck or an expansion) at time T in the past, which changed its size from N1 to N0 following an exponential model. The coalescent simulations were conditioned with various prior distributions to assess the lack of dependency of the posterior estimates of the three inferred parameters (N0, N1 and T). The prior and hyperpriors used in Msvar were showed in Additional file 2: Table S2. The prior distributions tested consisted of (i) a model in which both N0 and N1 had the same prior upper bound, (ii) where N0 had a wider prior than N1 (i.e. an expansion was allowed) and (iii) where N1 had a wider prior than N0 (i.e. a bottleneck was allowed). In addition, various upper bound values for the prior of the time of change were also tested. As to the microsatellite mutation rate, we used an average vertebrate microsatellite mutation rate of 10− 4 [48] allowed to vary across our markers between approximately 10− 2 and 10− 6 (Additional file 2: Table S2). The generation time of 6 years was selected based on previous results [49]. Each Msvar run consisted of 2 × 109 iterations of the MCMC algorithm discarding the first 10% of the coalescent simulations as burn-in. Convergence of the chains was assessed with Gelman and Rubin’s diagnostic [50] using the CODA library [51] implemented in R based on the four runs performed for each populations with different prior distributions.

Niche reconstruction

We use ENM to simulate the distribution pattern of P. shangchengensis under different climate backgrounds using climatic data and GPS coordinates of 30 distribution localities using the maximum entropy algorithm in the program MAXENT v. 3.3.3 k [52]. For the current climate predictions, we downloaded raster coverages of 19 environmental climatic variables from the WorldClim database (http://www.worldclim.org) at 30 arc-seconds resolution (~ 1 km2) and clipped these coverages to a region that encompassed the entire Dabie Mountains’s range (30 ~ 32°N, 114 ~ 117°E) [53]. We performed a correlation analysis between the nineteen climatic variables with ArcGIS 10.0 (ESRI) in the defined area, and then selected a subset of eight variables with a correlation ® lower than 0.85: BIO3 (Isothermality), BIO5 (Max Temperature of Warmest Month), BIO8 (Mean Temperature of Wettest Quarter), BIO9 (Mean Temperature of Driest Quarter), BIO13 (Precipitation of Wettest Month), BIO14 (Precipitation of Driest Month), BIO15 (Precipitation Seasonality), BIO18 (Precipitation of Warmest Quarter). Assuming niche conservatism over time [54], we projected the species distribution model to the climatic data layers for current environmental conditions, and for layers describing climatic variables during the mid-Holocene (6000 years ago), the last glacial maximum (LGM, about 22,000 years ago) and the last inter-glacial (LIG, 0.13 Mya ~ 0.12 Mya) [55], to obtain the predicted species distribution across time. For the LIG, LGM and the mid-Holocene predictions we used the same set of bioclimatic variables at a maximum resolution (2.5 arc minutes) was collected from the WorldClim database. The community climate system model (CCSM) [56] were used to reconstruct the paleodistribution.

We performed 15 replicates for the model with subsamples, 10,000 background points and 25% of points set for random testing. Analyses were run using default program parameter cumulative outputs, a 0.00001 convergence threshold, and a maximum of 5000 iterations. The area under the curve (AUC) of the receiver operating characteristic (ROC) plot were used for model evaluation [52]. All models were post processed and visualized in ArcGIS 10.0 (Esri, Redlands, CA, USA).

Results

Genetic diversity

We sequenced 984 bp of the ND2 gene and 1127 bp of the Cyt b gene from 193 individuals. Total Hd and π across all populations were 0.990 and 0.0345, respectively (Table 1). The TTZ and MW populations exhibited the highest haplotype diversity (Hd = 0.983) while the KJY population possessed the lowest value (Hd = 0.746). The BYM population exhibited the highest nucleotide diversity (π = 0.0030) while KHJ presented the lowest one (π = 0.0015).

For microsatellite data, 345 individuals were successfully genotyped at twelve loci (Table 1). No consistent departure from HWE or linkage disequilibrium was observed across the six populations. The nuclear genetic diversity varied between populations and summary statistic used with Na ranging from 8.83 (JTX) to 17.42 (BYM), HO ranging from 0.596 (KHJ) to 0.756 (BYM), while HE varied from 0.699 (KHJ) to 0.837 (BYM). Overall, among the six populations, BYM showed higher genetic diversity than the other five populations (Table 1).

Phylogenetic reconstructions

Partitionfinder identified that the most suitable models of evolution for the two mitochondrial genes studied here were HKY + I model for the combination of the first and second codon positions of both genes, while the GTR + G model for the third codon positions. The phylogenetic tree was split into two major branches, which corresponded to the northern two populations (JTX, KHJ) and southern four populations (MW, TTZ, BYM, KJY) respectively (Fig. 2). In the northern branch, the individuals from JTX and KHJ populations formed two distinct clades respectively, which acted as sister groups for each other. The southern branch was composed of three clades, with the individuals of MW and TTZ populations forming distinct monophyletic clusters respectively. The third southern clade was formed by the individuals of the BYM and KJY populations, however, these could not be completely differentiated form each other as the clade with all the BYM individuals also included some KJY individuals (Fig. 2). The network analysis recapitulated the phylogenetic tree and had five haplogroups with no shared haplotypes among different populations except for BYM and KJY. These five haplogroups are separated by multiple mutations (Fig. 2). The estimated time to the most recent common ancestor (TMRCA) of the entire ingroup was 3.96 Mya (95% HPD, 2.56–5.51 Mya), with the five lineages originating between 2.93 and 1.4 Mya (Fig. 2).

Fig. 2
figure 2

Phylogenetic relationships and haplotype network based on the mtDNA haplotypes of P. shangchengensis. a Phylogenetic tree of P. shangchengensis based on mitochondrial data. The numbers above branches are Bayesian Posterior Probabilities (BPP) indicating node support and the numbers underneath branches are split time estimates (Mya) with their 95% highest posterior density in square brackets; b Median-joining network with node sizes proportional by area to the frequencies of haplotypes. The numbers of mutations separating the haplotypes are shown on the branches, except for the one-step mutations. The little red diamond nodes indicate undetected haplotypes

Full size image

Population structure

The FST values among populations were all significant and ranged from 0.494 (BYM and KJY) to 0.962 (JTX and MW) for mtDNA data and from 0.087 (MW and TTZ) to 0.590 (JTX and KJY) for the microsatellite data (Table 2), indicating moderate to high level of genetic differentiation among populations. STRUCTURE analyses revealed a maximum ΔK value for K = 2 supporting the split into two main branches observed with the mitochondrial DNA (Fig. 3). However, it is worth noting that the STRUCTURE results showed a maximum posterior probability for K = 3 (Fig. 3) supporting the split of the southern branch observed with the mitochondrial data into a cluster formed by the MW and TTZ individuals, and another one formed by the BYM and KJY individuals. These two southern clusters share a substantial amount of the genetic variation as reflected by the admixture between individuals in both clusters (Fig. 3). Further splitting of the samples with STRUCTURE for K values of 4 to 6 does not result in better clustering results (Additional file 3: Figure S1). The inclusion of the geographical prior into the spatial model of Geneland with the microsatellite data consistently identified six groups in the ten replicas of the analysis carried out and which correspond to the six populations in this study (Fig. 4 and Additional file 4: Table S4).

Table 2 Pairwise FST among six populations of P. shangchengensis based on mtDNA sequence (below diagonal) and microsatellite data (above diagonal), respectively
Full size table
Fig. 3
figure 3

STRUCTURE clustering results based on microsatellite genotype data of six P. shangchengensis populations. a The linear relationship between LnP(D) and the number of clusters (K); b ΔK values as a function of K based on 5 runs; c STRUCTURE output for K = 2 and 3

Full size image
Fig. 4
figure 4

Spatial distribution of each group defined by Geneland for K = 6. a Cluster 1, b Cluster 2, c Cluster 3, d Cluster 4, e Cluster 5, and f Cluster 6. Clusters are indicated by areas with different intensities of color. Lighter-colored areas indicate a higher probability that individuals belong to that cluster

Full size image

The AMOVA analysis of the mitochondrial data with the predefined grouping of the samples according to the three clusters suggested by STRUCTURE (i.e. JTX and KHJ, MW and TTZ, BYM and KJY) shows that 61.02% of the genetic variation occurs among groups, 33.83% among populations within groups and 5.15% within populations (Additional file 5: Table S3). The same analysis performed on the microsatellite data shows that 7.75% of the genetic variation occurs among groups, 2.8% among populations within groups, and 89.46% within populations (Additional file 5: Table S3).

Historical demography

The neutrality test (Tajima’s D and Fu’s Fs), mismatch distribution analysis and BSPs were conducted on the mtDNA sequence data and excluding the two KJY haplotypes clustered within the BYM population. Fu’s Fs and Tajima’s D were negative in all six populations, however, while Fs was significant for all populations except KJY, only MW and BYM presented a significant D (Table 1). This results are further supported by the unimodal mismatch distributions observed in four out of six populations (TTZ and KJY showed multimodal or nearly unimodal distributions; Additional file 6: Figure S2 and Additional file 7: Table S5), and by the BSP results that showed that all populations experienced a demographic expansion (Additional file 8: Figure S3).

The Msvar analysis based on microsatellite genotype data was run four times independently under various demographic models. The Gelman and Rubin statistic inferred from these runs were smaller than the upper threshold of 1.2 for all populations indicating convergence of the MCMC. The Msvar results for the southern populations supports a stable demography in these populations with N0 and N1 posterior distributions with their modes almost completely overlapping each other. Contrastingly, the northern populations (JTX and KHJ) have an N0 (mode ~ 4870 and 5959, respectively) almost 15 times smaller than their N1 (mode ~ 74,989 and 81,752, respectively) indicating that these populations passed through a drastic bottleneck approximately 6000 years ago (Fig. 5, Table 3).

Fig. 5
figure 5

Posterior distributions of the demographic parameters inferred with Msvar of six populations based on microsatellite data (a–f). Estimated posterior distributions of current (N0, blue curve) and ancestral (N1, red curve) effective population sizes and time since population change (T, black curve) on a logarithmic scale based on a generation time of 6 years

Full size image
Table 3 Parameters Msvar analysis of each population in P. shangchengensis. Posterior probabilities and 95% highest posterior intervals for the parameters inferred in Msvar
Full size table

Ecological niche modeling

The output of MaxEnt consists of a grid map with each cell having an index of suitability between 0 and 1. The predicted distribution of P. shangchengensis based on ecological niche modeling closely matched its actual distribution (Fig. 6d), differing significantly from random (AUC = 0.998; 25% of data points set for random testing). Most sampling localities in this study fell within the area of highest predicted suitability. Of the 8 environmental variables used to construct the niche model, BIO5 (Max Temperature of Warmest Month), BIO9 (Mean Temperature of Driest Quarter) and BIO8 (Mean Temperature of Wettest Quarter) had the highest contributions to the model (69.1, 18.3 and 5.5%, respectively) (Additional file 9: Table S6). During LIG (0.13 Mya ~ 0.12 Mya), the potential distribution was the smaller than that of other three periods (Fig. 6a), and the distribution during mid Holocene (about 6000 years before the present) (Fig. 6c) of P. shangchengensis was smaller than that during the LGM (about 22,000 years before) (Fig. 6b). Based on the predictions of ENMs, although the distribution area has changed under different historical periods and climatic conditions, the results showed long-term stability in the availability of suitable habitat and long-term isolation among population distribution area by unsuitable environment at low altitude in P. shangchengensis (Fig. 6a, b, c, d).

Fig. 6
figure 6

Ecological niche modeling of P. shangchengensis under (a) Last Inter-glacial conditions, (b) Last Glacial Maximum conditions, (c) Mid Holocene, (d) Current conditions. Location data used for modeling are indicated as yellow dots

Full size image

Discussion

Genetic diversity

Small and isolated populations for prolonged periods often encounter a number of genetic risks, such as the loss of genetic variation due to genetic drift, inbreeding the potential accumulation of weakly deleterious mutations [57,58,59]. In this study, we described P. shangchengensis populations that have been isolated from each other for more than one million years, however, each populations still harbored relatively high levels of genetic diversity when compared to other narrowly distributed salamanders, such as Batrachuperus tibetanus (Hynobidae, Caudata) [60] and various Plethodontids (Plethodontidae, Caudata) such as P. ouachitae, P. fourchensis and P. caddoensis [10, 13, 14], or even some widely distributed amphibians such as Hensel’s swamp frog (Pseudopaludicola falcipes) in South America’s Pampas or Southwest China’s Leishan spiny toad (Leptobrachium leishanense) [61,62,63]. Several theories have been put forward to explain the high genetic diversity pattern observed in some species, such as species’ historic levels of genetic diversity, larger extant population size, longer generation times, random mating, gene flow, no inbreeding and lack of bottlenecks [64]. In general, during the Quaternary, dramatic climate oscillation led to expansions and contractions of glaciated areas, having an important impact on the montane habitats and leading some of them to become hotspots of genetic diversity [11, 12]. Based on the ENM prediction (Fig. 6), the Dabie Mountains could have provided a stable habitat for P. shangchengensis under different climate condition during the LIG, LGM and mid-Holocene periods, probably through the presence of favorable environments in long-term and stable montane stream habitats that prevented the collapse of most of its populations and thus favored the maintenance of their genetic variation.

In this study, a relatively low genetic diversity was observed in the JTX and KHJ populations. The amount of genetic variation may be correlated with the suitable habitat size for a population [4, 12, 17, 18, 65,66,67], with JTX and KHJ showing smaller current and past suitable habitat sizes in comparison to the other southern populations of P. shangchengensis. However, this observation is probably exacerbated by the demographic history of these two populations that passed through a relatively recent and strong bottleneck during the mid-Holocene further leading to a decrease in their genetic variation in comparison to the southern populations (Fig. 5, Table 3).

Geographic structure

Species restricted to montane sky island habitats often show high levels of inter-population genetic divergence and unique patterns of genetic structure [14, 19, 68,69,70], a pattern magnified in species distributed in mid to high-latitude mountain regions [11, 12]. In such habitats at high elevation, populations are often separated from each other by unsuitable habitat at lower elevation areas due to niche conservatism [10, 13, 14, 16, 17, 71]. This separation was often facilitated by climatic contrasts at different altitudes or orogenesis with the lower elevation areas becoming geographic barriers to gene flow leading to habitat fragmentation [16]. Consequently, species that moved along the altitudinal gradient, probably as a consequence of climate changes, may have originally formed a continuous population at lower altitudes that becomes subdivided when entering the sky islands and over time results in monophyletic groups through genetic drift [16, 17]. For example, in western Arkansas (USA), the unique physiographic features of the Ouachita Mountains area, coupled with species response to climatic factors, drove lineage divergence in three Plethodon species (P. ouachitae, P. fourchensis and P. caddoensis) resulting in the classic phylogeographic structure associated with stream drainages and mountains [10, 13, 14]. B. tibetanus (Hynobiidae, Caudata) is a typical stream salamander in the Qinling Mountains showing a split into two deeply divergent lineages that formed in the past 3 to 4 Mya as a response to the orogenesis of Qinling Mountains during the late Cenozoic [60].

In the Quaternary, sharp climatic oscillations caused species at mid-latitudes to experience the most severe population expansions and contractions [6]. In the mid-latitude mountain areas, the changes were also reflected in altitude [11, 12]. However, in East China, the coastal mountains intercept moisture and heat, transported by monsoons from the ocean and provide relative climatic and ecological stability [72,73,74,75]. The Dabie Mountains, composed of a chain of ancient isolated low-middle elevation massifs (Fig. 1) and located in a mid-latitude area of East Asia, are believed to have been able to maintain a relatively stable climate over the last several million years [73, 76]. Based on the ENMs prediction results, although the distribution area of P. shangchengensis has changed over different historical periods and climatic conditions, there has always been extensive suitable habitat for different populations (Fig. 6), even during the LGM that marked the most drastic climate oscillations of the whole Pleistocene [77]. The ENM results show that in lower elevation areas, unsuitable habitat for P. shangchengensis always existed, acting as a strict and effective isolation barrier for this species with strict niche conservatism. Once these discontinuous sky islands were formed, deep inter-population genetic divergences of P. shangchengensis accumulated gradually. Consequently, it is likely that an ancient P. shangchengensis population became fragmented after entering the current species distribution range with unsuitable habitat characterizing the geographic space between sky islands and having played a major role in the divergence of these populations.

Population demography

Generally, mid-latitude regions, and especially mountainous areas, experienced the most significant climatic changes due to the fluvial and glacial-dominated conditions during the Quaternary [11, 12]. In the mid-latitude regions of Europe and North America, the distributions of organisms were affected by drastic climate oscillations and glacial cycles during this period [6, 78]. Compared with the climate in Europe and North America, the mid-latitude regions in East Asia area were a mosaic of mountains lower than 2000 m characterized by a relatively mild Pleistocene climate due to a lack of ice sheets in many areas [73, 75, 79]. Although previous studies have suggested the glaciations during the Quaternary in East Asia were responsible for the cyclic population expansions or contractions for several taxa [80,81,82], evidence is available that for others (e.g. spiny frogs – Quasipaa boulengeri [83]) these events had a little or no limited, if any, effect. In this study, based on P. shangchengensis mtDNA data a long-term demographic expansion since the LIG was observed (0.13 Mya ~ 0.12 Mya) (Additional file 8: Figure S3) [84]. This observation agrees with the ENM prediction models for the species that show a relatively smaller distribution area at the beginning of LIG and relatively wide range of suitable habitats during LGM (Fig. 6). The onset of the demographic expansion during the LIG and population stabilization during the last glaciation of P. shangchengensis are in accordance with the observations for many other species in East Asia, such as the tufted deer (Elaphodus cephalophus) [82], the Chinese Hwamei (Leucodioptron canorum) [85] and the Chinese bamboo partridge (B. thoracica) [81]. In mid-latitude Europe, various reptiles also experienced similar demographic changes, e.g. the Italian and Maltese wall lizards (Podarcis siculus and P. filfolensis, respectively) [77, 86]. Obviously, some temperate species in mid-latitude regions, especially those insensitive to environmental changes, underwent attenuated glacial population contractions, or may have passed through population expansion [77, 81, 82, 85, 86].

In the present study, all P. shangchengensis populations, except the two northern populations (JTX and KHJ), remained stable during the Holocene (Fig. 5). The two northern populations showed evidence of a strong population decline approximately 6000 years ago (Fig. 5, Table 3). The decline of those two populations may be attributed to drastic climatic changes in the mid-Holocene that affected the northern boundary of the Dabie mountains more strongly that the rest of this southern range. Several recent studies highlight the mid-Holocene as a period of particularly profound climate variation, expressed as the decline of land air temperatures across much of the globe [87,88,89]. Additionally, compared with the habitat of four southern populations (MW, TTZ, BYM, KJY), the KHJ and JTX populations occupy a relatively smaller effective habitat area (Fig. 1) making them more sensitive to environmental changes in their geographic range as smaller populations have higher risk of population decline or extinction [90, 91].

Conclusion

Mid-latitude mountains are unique environments where species may exhibit amazing genetic divergence and population demography even in narrow areas. In this study, we show how the genetic variation of P. shangchengensis has been shaped by a combination of factors including the region’s topography, historical climate conditions and micro-climatic zones. The niche conservatism of P. shangchengensis within the context of its distribution range and evolutionary history support the sky island effect in the Eastern China Dabie Mountains. Although the current genetic variation in the species does not show signs of a population decline in modern times, the fragility and susceptibility to modern climate change of the environment it inhabits argues in favor of developing a management strategy for this and other species experiencing the sky island effect in the region.

Abbreviations

BSP:

Bayesian skyline plots

BYM:

Baimajian-Yaoluoping-Mingtangshan

CSM:

Community Climate System Model

Cyt b :

Cytochrome b

ENM:

Ecological niche modelling

JTX:

Jiaoyuan-Tanghui-Xiaolongtan

KHJ:

Kangwangzhai-Huangbaishan- Jiufengjian

KJY:

Kujingyuan

LGM:

Last Glacial Maximum

LIG:

Last Inter-glacial

MCMC:

Markov Chain Monte Carlo

MW:

Mazongling-Wochuan

ND2 :

NADH dehydrogenase 2

TMRCA:

Time to most recent common ancestor

TTZ:

Tiantangzhai

References

  1. Benton MJ. The red queen and the court jester: species diversity and the role of biotic and abiotic factors through time. Science. 2009;323(5915):728–32.

    Article  CAS  PubMed  Google Scholar 

  2. Baker PA, Fritz SC, Dick CW, Eckert AJ, Horton BK, Manzoni S, Ribas CC, Garzione CN, Battisti DS. The emerging field of geogenomics: constraining geological problems with genetic data. Earth-Sci Rev. 2014;135:38–47.

    Article  Google Scholar 

  3. Gaston KJ. Global patterns in biodiversity. Nature. 2000;405:220–7.

    Article  CAS  PubMed  Google Scholar 

  4. Favre A, Päckert M, Pauls SU, Jähnig SC, Uhl D, Michalak I, Muellner-Riehl AN. The role of the uplift of the Qinghai-Tibetan plateau for the evolution of Tibetan biotas. Biol Rev. 2015;90(1):236–53.

    Article  PubMed  Google Scholar 

  5. Yang SJ, Dong HL, Lei FM. Phylogeography of regional fauna on the Tibetan plateau: a review. Prog Nat Sci. 2009;19(7):789–99.

    Article  CAS  Google Scholar 

  6. Hewitt GM. Genetic consequences of climatic oscillations in the quaternary. Philos T Royal Soc B. 2004;359(1442):183–95.

    Article  CAS  Google Scholar 

  7. Ikeda DH, Max TL, Allan GJ, Lau MK, Shuster SM, Whitham TG. Genetically informed ecological niche models improve climate change predictions. Glob Chang Biol. 2017;23(1):164–76.

    Article  PubMed  Google Scholar 

  8. Zhou SZ, Wang J, Xu LB, Wang XL, Colgan PM, Mickelson DM. Glacial advances in southeastern Tibet during late Quaternary and their implications for climatic changes. Quatern Int. 2010;218(1):58–66.

    Google Scholar 

  9. Ni G, Kern E, Dong YW, Li Q, Park JK. More than meets the eye: the barrier effect of the Yangtze River outflow. Mol Ecol. 2017;26:4591–602.

    Article  PubMed  Google Scholar 

  10. Shepard DB, Burbrink FT. Local-scale environmental variation generates highly divergent lineages associated with stream drainages in a terrestrial salamander, Plethodon caddoensis. Mol Phylogenet Evol. 2011;59(2):399–411.

    Article  PubMed  Google Scholar 

  11. Herman F, Champagnac JD. Plio-Pleistocene increase of erosion rates in mountain belts in response to climate change. Terra Nova. 2016;28(1):2–10.

    Article  Google Scholar 

  12. Beniston M. Environmental change in mountains and uplands. New York: Routledge; 2016.

    Book  Google Scholar 

  13. Shepard DB, Burbrink FT. Lineage diversification and historical demography of a sky island salamander, Plethodon ouachitae, from the interior highlands. Mol Ecol. 2008;17(24):5315–35.

    Article  PubMed  Google Scholar 

  14. Shepard DB, Burbrink FT. Phylogeographic and demographic effects of Pleistocene climatic fluctuations in a montane salamander, Plethodon fourchensis. Mol Ecol. 2009;18(10):2243–62.

    Article  PubMed  Google Scholar 

  15. Wiens JJ, Harrison R. Speciation and ecology revisited: phylogenetic niche conservatism and the origin of species. Evolution. 2004;58(1):193–7.

    Article  PubMed  Google Scholar 

  16. McCormack JE, Huang H, Knowles LL, Gillespie R, Clague D. Sky islands. Encyclopedia of Islands. 2009;4:841–3.

    Google Scholar 

  17. He K, Jiang XL. Sky islands of Southwest China. I: an overview of phylogeographic patterns. Chin Sci Bull. 2014;59(7):585–97.

    Article  Google Scholar 

  18. Taubmann J, Theissinger K, Feldheim KA, Laube I, Graf W, Haase P, Johannesen J, Pauls SU. Modelling range shifts and assessing genetic diversity distribution of the montane aquatic mayfly Ameletus inopinatus in Europe under climate change scenarios. Conserv Genet. 2011;12(2):503–15.

    Article  Google Scholar 

  19. Wu YK, Wang YZ, Jiang K, Hanken J. Significance of pre–quaternary climate change for montane species diversity: insights from Asian salamanders (Salamandridae: Pachytriton). Mol Phylogenet Evol. 2013;66(1):380–90.

    Article  PubMed  Google Scholar 

  20. Zheng YH, Zhang Y, Shao XM, Yin ZY, Zhang J. Temperature variability inferred from tree-ring widths in the Dabie Mountains of subtropical Central China. Trees. 2012;26(6):1887–94.

    Article  Google Scholar 

  21. Fei L, Ye CY, Jiang JP. Colored atlas of Chinese amphibians and their distributions. Chengdu: Sichuan Publishing House of Science and Technology; 2012.

    Google Scholar 

  22. Sambrook J, Fritsch EF, Maniatis T. Molecular cloning. New York: Cold Spring Harbor Laboratory Press; 1989.

    Google Scholar 

  23. Clarke KR, Gorley RN. Primer v5. Primer–E Ltd. 2001.

  24. Wang H, Shi WB, Zhu LF, Luo X, Chen Z, Xie WL, Han DM, Chang Q, Zhang BW. Isolation and characterization of fourteen polymorphic tetranucleotide microsatellites for Shangcheng stout salamander (Pachyhynobius shangchengensis). Conserv Genet Res. 2010;2(1):121–3.

    Article  Google Scholar 

  25. Holland MM, Parson W. GeneMarker® HID: a reliable software tool for the analysis of forensic STR data. J Forensic Sci. 2011;56(1):29–35.

    Article  PubMed  Google Scholar 

  26. Burland TG. DNASTAR’s Lasergene sequence analysis software. In: Misener S, Krawetz SA, editors. Bioinformatics methods and protocols. Totowa, NJ: Methods in Molecular Biology™; 2000. p. 71–91.

    Google Scholar 

  27. Larkin MA, Blackshields G, Brown NP, Chenna R, McGettigan PA, McWilliam H, Valentin F, Wallace IM, Wilm A, Lopez R, et al. Clustal W and Clustal X version 2.0. Bioinformatics. 2007;23(21):2947–8.

    Article  CAS  PubMed  Google Scholar 

  28. Van Oosterhout C, Hutchinson WF, Wills DPM, Shipley P. MICRO–CHECKER: software for identifying and correcting genotyping errors in microsatellite data. Mol Ecol Notes. 2004;4:535–8.

    Article  CAS  Google Scholar 

  29. Rousset F. Genepop’007: a complete re-implementation of the genepop software for windows and Linux. Mol Ecol Res. 2008;8(1):103–6.

    Article  Google Scholar 

  30. Peakall R, Smouse PE. GENALEX 6: genetic analysis in excel. Population genetic software for teaching and research. Mol Ecol Notes. 2006;6:288–95.

    Article  Google Scholar 

  31. Belkhir K, Borsa P, Chikhi L. 1996–2001 GENETIX 4.02, logiciel sous Windows TM pour la génétique des populations. 2001.

  32. Drummond AJ, Suchard MA, Xie D, Rambaut A. Bayesian phylogenetics with BEAUti and the BEAST 1.7. Mol Biol Evol. 2012;29(8):1969–73.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Lanfear R, Calcott B, Ho SYW, Guindon S. PartitionFinder: combined selection of partitioning schemes and substitution models for phylogenetic analyses. Mol Biol Evol. 2012;29:1695–701.

    Article  CAS  PubMed  Google Scholar 

  34. Chen MY, Mao RL, Liang D, Kuro-O M, Zeng XM, Zhang P. A reinvestigation of phylogeny and divergence times of Hynobiidae (Amphibia, Caudata) based on 29 nuclear genes. Mol Phylogenet Evol. 2015;83:1–6.

    Article  PubMed  Google Scholar 

  35. Rambaut A, Suchard M, Xie D, Drummond AJ. MCMC trace analysis package (version 1.6); 2014. Available at: http://tree.bio.ed.ac.uk/software/tracer/

  36. Rambaut A, Drummond AJ. TreeAnnotator version 1.8; 2016. Available at: http://beast.community/programs

  37. Rambaut A. FigTree v1.4.3; 2016. Available at: http://tree.bio.ed.ac.uk/software/figtree/

  38. Bandelt HJ, Forster P, Röhl A. Median–joining networks for inferring intraspecific phylogenies. Mol Biol Evol. 1999;16(1):37–48.

    Article  CAS  PubMed  Google Scholar 

  39. Excoffier L, Laval G, Schneider S. Arlequinver 3.0: an integrated software package for population genetic analysis. Evol Bioinforma. 2005;1:47–50.

    Article  CAS  Google Scholar 

  40. Pritchard JK, Stephens M, Donnelly P. Inference of population structure using multilocus genotype data. Genetics. 2000;155(2):945–59.

    CAS  PubMed  PubMed Central  Google Scholar 

  41. Evanno G, Regnaut S, Goudet J. Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Mol Ecol. 2005;14(8):2611–20.

    Article  CAS  PubMed  Google Scholar 

  42. Rosenberg NA. DISTRUCT: a program for the graphical display of population structure. Mol Ecol Notes. 2004;4(1):137–8.

    Article  Google Scholar 

  43. Guillot G, Mortier F, Estoup A. Geneland: a computer package for landscape genetics. Mol Ecol Resour. 2005;5(3):712–5.

    Article  CAS  Google Scholar 

  44. Harpending HC. Signature of ancient population growth in a low-resolution mitochondrial DNA mismatch distribution. Hum Biol. 1994:591–600.

  45. Tajima F. The effect of change in population size on DNA polymorphism. Genetics. 1989;123(3):597–601.

    CAS  PubMed  PubMed Central  Google Scholar 

  46. Fu YX. Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics. 1997;147(2):915–25.

    CAS  PubMed  PubMed Central  Google Scholar 

  47. Storz JF, Beaumont MA. Testing for genetic evidence of population expansion and contraction: an empirical analysis of microsatellite DNA variation using a hierarchical Bayesian model. Evolution. 2002;56:154–66.

    Article  CAS  PubMed  Google Scholar 

  48. Bulut Z, McCormick CR, Gopurenko D, Williams RN, Bos DH, DeWoody JA. Microsatellite mutation rates in the eastern tiger salamander (Ambystoma tigrinum tigrinum) differ 10-fold across loci. Genetica. 2009;136(3):501–4.

    Article  CAS  PubMed  Google Scholar 

  49. Pasmans F, Janssens GPJ, Sparreboom M, Jiang JP, Nishikawa K. Reproduction, development, and growth response to captive diets in the Shangcheng stout salamander, Pachyhynobius shangchengensis (Amphibia, Urodela, Hynobiidae). Asian Herpetol Res. 2012;03(3):192–7.

    Article  Google Scholar 

  50. Brooks S, Gerlman A. General methods for monitoring convergence of iterative simulations. J Comput Graph Stat. 1998;7:434–55.

    Google Scholar 

  51. Plummer M, Best N, Cowles K, Vines K. Coda: output analysis and diagnostic for MCMC. R News. 2006;6:7–11.

    Google Scholar 

  52. Phillips SJ, Anderson RP, Schapire RE. Maximum entropy modeling of species geographic distributions. Ecol Model. 2006;190(3):231–59.

    Article  Google Scholar 

  53. Hijmans RJ, Cameron S, Parra J, Jones PG, Jarvis A. WorldClim, version 1.3. Berkeley: University of California; 2005.

    Google Scholar 

  54. Wiens JJ, Graham CH. Niche conservatism: integrating evolution, eology, and conservation biology. Annu Rev Ecol Evol S. 2005;36(1):519–39.

    Article  Google Scholar 

  55. Ottobliesner BL, Marshall SJ, Overpeck JT, Miller GH, Hu A. Simulating Arctic climate warmth and icefield retreat in the last interglaciation. Science. 2006;311(5768):1751–3.

    Article  CAS  Google Scholar 

  56. Collins WD, Bitz CM, Blackmon ML, Bonan GB, Bretherton CS, Carton JA, Chang P, Doney SC, Hack JJ, Henderson TB. The community climate system model version 3 (CCSM3). B Am Meteorol Soc. 2003;82(11):2357–76.

    Google Scholar 

  57. Keller LF, Waller DM. Inbreeding effects in wild populations. Trends Ecol Evol. 2002;17(5):230–41.

    Article  Google Scholar 

  58. Lynch M, Conery J, Burger R. Mutation accumulation and the extinction of small populations. Am Nat. 1995;146(4):489–518.

    Article  Google Scholar 

  59. Briskie JV, Mackintosh M. Hatching failure increases with severity of population bottlenecks in birds. Proc Natl Acad Sci U S A. 2004;101(2):558–61.

    Article  CAS  PubMed  Google Scholar 

  60. Huang ZS, Yu FL, Gong HS, Song YL, Zeng ZG, Zhang Q. Phylogeographical structure and demographic expansion in the endemic alpine stream salamander (Hynobiidae: Batrachuperus) of the Qinling Mountains. Sci Rep. 2017;7(1871):1–13.

    Google Scholar 

  61. Wang E, Wijk REV, Braun MS, Wink M. Gene flow and genetic drift contribute to high genetic diversity with low phylogeographical structure in European hoopoes (Upupa epops). Mol Phylogenet Evol. 2017;113:113–25.

    Article  PubMed  Google Scholar 

  62. Langone JA, Camargo A, de Sá RO. High genetic diversity but low population structure in the frog Pseudopaludicola falcipes (Hensel, 1867) (Amphibia, Anura) from the pampas of South America. Mol Phylogenet Evol 2016; 95(4):137–151.

  63. Zhang W, Luo ZH, Zhao M, Wu H. High genetic diversity in the endangered and narrowly distributed amphibian species Leptobrachium leishanense. Integr Zool. 2015;10(5):465–81.

    Article  PubMed  Google Scholar 

  64. Allentoft ME, O’Brien J. Global amphibian declines, loss of genetic diversity and fitness: a review. Diversity. 2010;2(1):47–71.

    Article  Google Scholar 

  65. Villard M, Metzger JP. Beyond the fragmentation debate: a conceptual model to predict when habitat configuration really matters. J Appl Ecol. 2014;51(2):309–18.

    Article  Google Scholar 

  66. Petit RJ, Aguinagalde I, de Beaulieu JL, Bittkau C, Brewer S, Cheddadi R, Ennos R, Fineschi S, Grivet D, Lascoux M. Glacial refugia: hotspots but not melting pots of genetic diversity. Science. 2003;300(5625):1563–5.

    Article  CAS  PubMed  Google Scholar 

  67. Rull V. Microrefugia. J Biogeogr. 2009;36(3):481–4.

    Article  Google Scholar 

  68. Kozak KH, Wiens JJ. Does niche conservatism promote speciation? A case study in north American salamanders. Evolution. 2006;60(12):2604–21.

    Article  PubMed  Google Scholar 

  69. Pauls SU, Lumbsch T, Haase P. Phylogeography of the montane caddisfly Drusus discolor: evidence for multiple refugia and periglacial survival. Mol Ecol. 2006;15:2153–69.

    Article  CAS  PubMed  Google Scholar 

  70. Zhu LF, Zhang SN, Gu XD, Wei FW. Significant genetic boundaries and spatial dynamics of giant pandas occupying fragmented habitat across Southwest China. Mol Ecol. 2011;20(6):1122–32.

    Article  PubMed  Google Scholar 

  71. Cox SC, RP P–J, Habel JC, Amakobe BA, Day JJ. Niche divergence promotes rapid diversification of east African sky island white-eyes (Aves: Zosteropidae). Mol Ecol. 2014;23(16):4103–18.

    Article  PubMed  PubMed Central  Google Scholar 

  72. An ZS, Kutzbach JE, Prell WL, Porter SC. Evolution of Asian monsoons and phased uplift of the Himalaya–Tibetan plateau since Late Miocene times. Nature. 2001;411(6833):62–6.

    Article  CAS  Google Scholar 

  73. Ju L, Wang H, Jiang D. Simulation of the last glacial maximum climate over East Asia with a regional climate model nested in a general circulation model. Palaeogeogr Palaeocl. 2007;248:376–90.

    Article  Google Scholar 

  74. Sun XG, Wang PX. How old is the Asian monsoon system?—Palaeobotanical records from China. Palaeogeogra Palaeocl. 2005;222(3):181–222.

    Article  Google Scholar 

  75. Weaver AJ, Eby M, Fanning AF, Wiebe EC. Simulated influence of carbon dioxide, orbital forcing and ice sheets on the climate of the last glacial maximum. Nature. 1998;394(6696):847–53.

    Article  CAS  Google Scholar 

  76. Zhao B, Zhang J, Sun XH. Eco-environmental vulnerability evaluation based on GIS in Tongbai–Dabie Mountain area of Huai River basin. Res Soil Water Conserv. 2009;16(3):135–8.

    Google Scholar 

  77. Senczuk G, Colangelo P, Simone ED, Aloise G, Castiglia R. A combination of long term fragmentation and glacial persistence drove the evolutionary history of the Italian wall lizard Podarcis siculus. BMC Evol Biol. 2017;17(1):1–15.

    Article  CAS  Google Scholar 

  78. Schmitt T. Molecular biogeography of Europe: Pleistocene cycles and postglacial trends. Front Zool. 2007;4(1):1–13.

    Article  CAS  Google Scholar 

  79. Rost KT. Paleoclimatic field studies in and along the Qinling Shan (Central China). GeoJournal. 1994;34:107–20.

    Article  Google Scholar 

  80. Ding L, Gan XN, He SP, Zhao EM. 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. 2011;20(9):1905–22.

    Article  PubMed  Google Scholar 

  81. Huang ZH, Liu NF, Liang W, Zhang YY, Liao XJ, Ruan LZ, Yang ZS. Phylogeography of Chinese bamboo partridge, Bambusicola thoracica thoracica (Aves: Galliformes) in South China: inference from mitochondrial DNA control-region sequences. Mol Phylogenet Evol. 2010;56(1):273–80.

    Article  PubMed  Google Scholar 

  82. Sun ZL, Pan T, Wang H, Pang MJ, Zhang BW. Yangtze River, an insignificant genetic boundary in tufted deer (Elaphodus cephalophus): the evidence from a first population genetics study. PeerJ. 2016;4(9):e2654.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  83. Yan F, Zhou WW, Zhao HT, Yuan ZY, Wang YY, Jiang K, Jin JQ, Murphy RW, Che J, Zhang YP. Geological events play a larger role than Pleistocene climatic fluctuations in driving the genetic structure of Quasipaa boulengeri (Anura: Dicroglossidae). Mol Ecol. 2013;22(4):1120–33.

    Article  PubMed  Google Scholar 

  84. Yuan D, Cheng H, Edwards RL, Dykoski CA, Kelly MJ, Zhang M, Qing J, Lin Y, Wang Y, Wu J. Timing, duration, and transitions of the last interglacial Asian monsoon. Science. 2004;304(5670):575–8.

    Article  CAS  PubMed  Google Scholar 

  85. Li SH, Yeung CKL, Feinstein J, Han LX, Manh HL, Wang CX, Ding P. 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. 2009;18(4):622–33.

    Article  CAS  PubMed  Google Scholar 

  86. Salvi D, Schembri PJ, Sciberras A, Harris DJ. Evolutionary history of the Maltese wall lizard Podarcis filfolensis: insights on the ‘expansion–contraction’ model of Pleistocene biogeography. Mol Ecol. 2014;23(5):1167–87.

    Article  PubMed  Google Scholar 

  87. Cheng H, Fleitmann D, Edwards LR, Wang X, Cruz FW, Auler AS, Mangini A, Wang Y, Kong X, Burns SJ, et al. Timing and structure of the 8.2 kyr B.P. event inferred from 18O records of stalagmites from China, Oman, and Brazil. Geology. 2009;37(11):1007–10.

    Article  CAS  Google Scholar 

  88. Wanner H, Solomina ON, Grosjean M, Ritz SP, Jetel M. Structure and origin of Holocene cold events. Quaternary Sci Rev. 2011;30(21):3109–23.

    Article  Google Scholar 

  89. Bond G, Bonani G. Persistent solar influence on North Atlantic climate during the Holocene. Science. 2001;294(5549):2130–6.

    Article  CAS  PubMed  Google Scholar 

  90. Primack RB. A primer of conservation biology. Sunderland: Sinauer Associates; 1995.

    Google Scholar 

  91. Green DM. The ecology of extinction: population fluctuation and decline in amphibians. Biol Conserv. 2003;111(3):331–43.

    Article  Google Scholar 

Download references

Acknowledgements

We thank Dr. Wen-Liang Zhou and Dr. Kai Zhao for their help in sample collecting, Prof. Mike Bruford (Cardiff University) and Dr. Xing-Guang Li for data analyses. We thank Tianma National Nature Reserve and Yaoluoping Nature Reserve for their help in the Biodiversity Survey and sample collection. We thank Dr. Raul E. Diaz (Southeastern Louisiana University) for reviewing our manuscript and improving the English language. We thank three anonymous reviewer for their suggestions.

Funding

This work was supported by the National Natural Science Foundation of China (Grant No. 31272332, 30800116) and Anhui Province Higher Education Revitalization Plan (2014 Colleges and Universities Outstanding Youth Talent Support Program) in the design of study, field work, sample collection, lab work and data analysis. The writing and publishing was supported by Anhui Province Academic and Technical Leader & Backup Candidate Academic Research Activities Fund (2017H130).

Availability of data and materials

Sequences are deposited in GenBank (MF581931–MF582316).

Author information

Author notes

  1. Tao Pan, Hui Wang and Pablo Orozcoterwengel contributed equally to this work.

Authors and Affiliations

  1. Anhui Key Laboratory of Eco–engineering and Bio-technique, School of Life Sciences, Anhui University, Hefei, 230601, Anhui, China

    Tao Pan, Hui Wang, Gui-You Wu, Li-Fu Qian, Zhong-Lou Sun, Wen-Bo Shi & Bao-Wei Zhang

  2. School of Life Sciences, Anhui Normal University, Wuhu, 241000, Anhui, China

    Tao Pan, Peng Yan & Xiao-Bing Wu

  3. School of Biosciences, Cardiff University, Cardiff, UK

    Pablo Orozcoterwengel & Bao-Wei Zhang

  4. Analytical and Testing Center, Nanjing Normal University, Nanjing, 210046, Jiangsu, China

    Chao-Chao Hu

Authors
  1. Tao Pan
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Hui Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Pablo Orozcoterwengel
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Chao-Chao Hu
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Gui-You Wu
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Li-Fu Qian
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Zhong-Lou Sun
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Wen-Bo Shi
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Peng Yan
    View author publications

    You can also search for this author in PubMed Google Scholar

  10. Xiao-Bing Wu
    View author publications

    You can also search for this author in PubMed Google Scholar

  11. Bao-Wei Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

BWZ led the research team. BWZ, XBW, TP, HW and PO designed the research. TP, HW, LFQ, ZLS, WBS and PY carried out sample collection. TP, HW, LFQ, ZLS, WBS and PY performed lab work. TP, HW, PO, CCH, GYW and PY analyzed data. TP, HW and PO wrote the paper. All authors have read and approved the manuscript for submission.

Corresponding author

Correspondence to Bao-Wei Zhang.

Ethics declarations

Ethics approval and consent to participate

In the present study, our experimental procedures and sample collection complied with the current laws on animal welfare and research in China, and were specifically approved by the Animal Research Ethics Committee of Anhui Normal University.

Consent for publication

Not applicable.

Competing interests

The authors’ declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1:

Table S1. Details on the primers for amplified the ND2 gene, Cyt b gene and twelve microsatellite loci of P. shangchengensis. (DOCX 19 kb)

Additional file 2:

Table S2. Prior and hyperpriors used in Msvar based on microsatellite data. (DOCX 17 kb)

Additional file 3:

Figure S1. STRUCTURE clustering results for K = 4, 5 and 6. (TIF 2521 kb)

Additional file 4:

Table S4. Multiple runs for inferring the number of populations of P. shangchengensis with Geneland. Bold indicates the highest average posterior probability. (DOCX 16 kb)

Additional file 5:

Table S3. Results of hierarchical AMOVA based on mtDNA and microsatellites. Groups are set as JTX and KHJ, MW and TTZ, BYM and KJY. (DOCX 17 kb)

Additional file 6:

Figure S2. Mismatch distributions analyses for the six populations of P. shangchengensis based on the mtDNA sequence data (a–f). The two KJY individuals clustering with BYM were excluded from the demographic history analyses. The x coordinate represents the number of differences between pairs of sequences and the y coordinate represents the frequencies of pairwise differences. The blue histograms are the observed frequencies of pairwise divergences among sequences and the green line refers to the expected shape of the distribution under the model of population expansion. (TIF 2468 kb)

Additional file 7:

Table S5. Results of demographic statistics of P. shangchengensis based on mtDNA data. (DOCX 18 kb)

Additional file 8:

Figure S3. Bayesian skyline plot of effective population size in P. shangchengensis based on the mtDNA data (a–f). The BSP analyses of the KJY population was calculated excluding the individuals with two haplotypes clustered more closely to the BYM population. The x-axis indicates time in Mya BP, and the y-axis indicates the effective population size in units of Neτ (the product of effective population size and generation time in Mya). The blue areas represent 95% highest posterior density. Time is expressed in million years. (TIF 2011 kb)

Additional file 9:

Table S6. Results from Principal Components Analysis on environmental variables used in comparison of environments among occurrence localities for P. shangchengensis and the relative contribution of each variable to the niche model. (DOCX 16 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Pan, T., Wang, H., Orozcoterwengel, P. et al. Long-term sky islands generate highly divergent lineages of a narrowly distributed stream salamander (Pachyhynobius shangchengensis) in mid-latitude mountains of East Asia. BMC Evol Biol 19, 1 (2019). https://doi.org/10.1186/s12862-018-1333-8

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12862-018-1333-8

Keywords

  • Pachyhynobius shangchengensis
  • Sky island
  • Genetic diversity
  • Phylogeography
  • Population demography
  • Niche conservatism

BMC Ecology and Evolution

ISSN: 2730-7182

Contact us