Rapid genetic divergence and mitonuclear discordance in the Taliang knobby newt ( Liangshantriton taliangensis, Salamandridae, Caudata) and their driving forces.

The Hengduan Mountains Region (HMR) is the largest "evolutionary frontier" of the northern temperate zone, and the origin and maintenance of species in this area is a research hotspot. Exploring species-specific responses to historical and contemporary environmental changes will improve our understanding of the role of this region in maintaining biodiversity. In this study, mitochondrial and microsatellite diversities were used to assess the contributions of paleogeological events, Pleistocene climatic oscillations, and contemporary landscape characteristics to the rapid intraspecific diversification of Liangshantriton taliangensis, a vulnerable amphibian species endemic to several sky-island mountains in the southeastern HMR. Divergence date estimations suggested that the East Asian monsoon, local uplifting events (Xigeda Formation strata), and Early-Middle Pleistocene transition (EMPT) promoted rapid divergence of L. taliangensis during the Pleistocene, yielding eight mitochondrial lineages and six nuclear genetic lineages. Moreover, population genetic structures were mainly fixed through isolation by resistance. Multiple in situ refugia were identified by ecological niche models and high genetic diversity, which played crucial roles in the persistence and divergence of L. taliangensis during glacial-interglacial cycles. Dramatic climatic fluctuations further promoted recurrent isolation and admixing of populations in scattered glacial refugia. The apparent mitonuclear discordance was likely the result of introgression by secondary contact and/or female-biased dispersal. Postglacial expansion generated two major secondary contact zones (Ganluo (GL) and Chuhongjue (CHJ)). Identification of conservation management units and dispersal corridors offers important recommendations for the conservation of this species.

INTRODUCTION

Pleistocene climatic oscillations and Meso-Cenozoic intracontinental orogenesis are considered essential drivers of the current geographical distribution and genetic divergence of organisms (Favre et al., 2015; Gillespie & Roderick, 2014; Rana et al., 2021). Topographical changes, such as mountain and river systems caused by orogenesis, are commonly regarded as geographical barriers that promote biodiversity, especially in species with limited dispersal abilities (Rana et al., 2021; Sánchez-Montes et al., 2018). Mountainous regions often become hotspots of species diversity and lineage divergence (Kozak & Wiens, 2010; Noroozi et al., 2018), with the uplift of the Qinghai-Xizang (Tibet) Plateau (QTP) receiving considerable attention (He et al., 2021). The barrier effect caused by this uplift has been documented for various fauna and flora taxa from the QTP and adjacent areas (Huang et al., 2017; Wang et al., 2017; Ye et al., 2019). However, the geological history of these regions is highly complex, and many details on the timing and patterns of geological uplift remain controversial (Harrison et al., 1992; Li & Fang, 1999; Renner, 2016; Spicer et al., 2021; Su et al., 2019). Thus, interpretation of the effect of QTP uplift on species formation and lineage divergence needs to focus on differences in biota from the QTP and surrounding regions (Xing & Ree, 2017). This is particularly important for relatively young species that developed from the Miocene onward. Uplift of the QTP during the Miocene further contributed to the establishment of the South Asian and East Asian monsoon climate (An et al., 2001; Liu & Dong, 2013). Furthermore, QTP uplift shaped the warm and humid environment, dominated by tropical and subtropical vegetation on the southeastern edge, especially the Hengduan Mountains Region (HMR). Many unoccupied ecological niches formed by new habitats in the region promoted the evolutionary radiation of organisms and enhanced their diversification potential (Favre et al., 2015; Wang et al., 2015, 2018; Xing & Ree, 2017).

Climatic fluctuations in the Quaternary caused repeated and extreme environmental changes, which profoundly affected the distribution, dispersal, and adaptation of organisms (Davis & Shaw, 2001; Hewitt, 2000). Different distributional and physiological characteristics allowed species to respond to glaciation differently (Huang et al., 2017; Qu et al., 2010; Yu et al., 2013). Alpine species adapted to temperate and warm habitats retreated southward to low-altitude refugia during glacial periods, thus avoiding glacial cover and the cold-dry climate, and expanded and recolonized the northern and plateau areas when the climate again warmed (Gutiérrez-Rodríguez et al., 2017; Qiu et al., 2011). Cold-tolerant species or species with distribution areas in climatically stable zones often survivedin situ instead of migrating to refugia (Sánchez-Montes et al., 2019; Wang et al., 2017). Microrefugia play an important role in the rapid divergence and persistence of lineages in mountainous regions with severe climatic fluctuations (He et al., 2021; Rana et al., 2021; Vasconcellos et al., 2019). Furthermore, shifts in species range caused by glacial-interglacial cycles provide recurrent opportunities for expansion and admixture of lineages, thus facilitating the occurrence of mitonuclear discordance (Dufresnes et al., 2020; Phuong et al., 2017).

The HMR is a transition region between the eastern circum-Pacific belt and the western ancient Mediterranean belt and is a constitutional part of the QTP. The area (~ N23°–33° and E94°–103°) is characterized by a series of parallel high mountains and great rivers running north to south. The terrain becomes gentle from the northwest to southeast, but with sharp altitudinal differences. The HMR, Yunnan-Guizhou Plateau, and Bashan Mountains are known as “the southwest China sky-island complex” (He & Jiang, 2014). Dispersal ability and habitat preference may affect the response patterns of species to landscape features (Robertson et al., 2018). The mountains and rivers of the sky-island complex constitute a potential barrier for gene flow and thus may promote rapid species diversification (Atlas & Fu, 2019; Deng et al., 2020; Pan et al., 2020) or may provide dispersal corridors for other species (Qu et al., 2011; Zhang et al., 2010). The southeastern fringe of the HMR consists of large ice-free areas, especially at low altitudes, which are considered as glacial refugia for a variety of plants and animals (Xie et al., 2019; Zhan et al., 2011). The HMR has a unique geographical history (Xing & Ree, 2017). Notably, while the area likely uplifted during the late Miocene to Pliocene (Favre et al., 2015), recent fossil evidence indicates that part of the HMR may have uplifted earlier, i.e., during the early Oligocene (Su et al., 2019). Evolutionary radiation of sky-island species can often be traced to the Pleistocene or earlier, especially for endemic species (He & Jiang, 2014; Yu et al., 2013). The HMR plays an important role in maintaining biodiversity and is considered the largest “evolutionary frontier” of the northern temperate zone (Boufford, 2014; López‐Pujol et al., 2011) and one of the most vulnerable biodiversity hotspots due to climate change and human disturbance (Ding et al., 2020; Myers et al., 2000). The inter-relationships between specific sky-island species responses and topography as well as climatic fluctuations in the HMR need to be further explored (He & Jiang, 2014). This is particularly important for amphibians with generally weak dispersal ability and more apparent genetic structure.

The Taliang knobby newt, Liangshantriton taliangensis Liu, 1950 belongs to Salamandridae Goldfuss, 1820 (Fei et al., 2012). The species is endemic to several sky-island mountains in the southeastern HMR and is distributed at altitudes from 1 390 to 3 000 m a.s.l. (Fei et al., 2012), where it usually inhabits ponds and wetlands with high humidity and heavy vegetation. Its distribution area is in the transition zone between the first and second steps of China’s terrain, containing a diverse topography and belonging to the branches of the Daxueshan Mountains. The L. taliangensis population continues to decline and is threatened by habitat degradation and excessive collection. Consequently, it has been classified as vulnerable by the International Union for Conservation of Nature red list (IUCN SSC Amphibian Specialist Group, 2020) and the red list of China’s vertebrates (Jiang et al., 2016). Research on the origin and lineage divergence of L. taliangensis not only improves our understanding of the interaction mechanism between historical and contemporary environmental factors and sky-island species divergence in the southwest mountainous area in China (especially the HMR), but also provides a reference for the conservation of L. taliangensis.

In this study, both mitochondrial and microsatellite markers were used to assess the rapid lineage divergence and genetic structure of L. taliangensis populations in southeastern HMR. The ecological niche model (ENM) and landscape genetics (isolation-by-resistance (IBR) and isolation-by-environment (IBE)) were used to evaluate the role of regional uplift events, Quaternary climatic oscillations, and landscape factors on intraspecific diversification. Based on the study objectives, we: (a) explored the phylogenetic relationship and genetic structure of L. taliangensis; (b) evaluated the historical and contemporary landscape and climatic factors that contributed to the rapid divergence of L. taliangensis; (c) identified the responses of L. taliangensis to Pleistocene glaciation; and (d) suggested conservation implications for L. taliangensis.

MATERIALS AND METHODS

Sampling, DNA extraction, and amplification

Tissue samples (426) of L. taliangensis were collected from 16 sites throughout the mountain ranges of the southeastern HMR in China, including the Gonggar (GG), Xiaoxiangling (XXL), and Liangshan mountains (LS), which include Xiaoliangshan (XLS) and Daliangshan (DLS) (Figure 1; Supplementary Table S1). All procedures followed the guidelines established by the Animal Care and Use Committee of the Chengdu Institute of Biology, Chinese Academy of Sciences (Permit No.: 2017-AR-JJP-03). Three mitochondrial genes were amplified, including NADH dehydrogenase subunit 2 (ND2), cytochrome b (cyt b), and cytochrome c oxidase subunit I (COI). Multiple primer pairs were used for each gene amplification (Supplementary Table S2). Nine species from four clades of Tylototriton sensu lato (s.l.) were used as outgroups to root the trees (Wang et al., 2018). Among these species, the DNA of Tylototriton pseudoverrucosus Hou, Gu, Zhang, Zeng, and Lu, 2012 was extracted from a specimen (Specimen No.: CIB100802), while all others were obtained from the GenBank database (Supplementary Table S3). Eleven high polymorphism microsatellite markers (DLY32,DLY39, DLY42,DLY74, DLY76, DLY96, DLY124, DLY131, DLY134, DLY148, and DLY153) were developed and used to characterize the genotypes of each individual (Chen et al., 2019; Shu, 2020). Further details are presented in the Supplementary Methods.

Figure 1

Results of the Bayesian clustering analysis across Liangshantriton taliangensis distribution range

Black dashed line represents four main distribution areas of L. taliangensis (Gonggar (GG), Xiaoxiangling (XXL), Xiaoliangshan (XLS), and Daliangshan mountains (DLS)). A: Geographic distribution of two main genetic clusters, as identified by microsatellite data. West cluster (red) contains populations of XXL and GG mountains (XXL-GG); East cluster (blue) contains populations of Liangshan Mountains (LS), including DLS and XLS. Black arrows indicate gene flow. B: Geographic distribution of four subclusters of West cluster. Subcluster A (orange) contains populations of GG mountains; subclusters B (yellow), C (green), and D (purple) belong to XXL mountains. C: Geographic distribution of two subclusters in East cluster. Subcluster E (magenta) contains DLS and GL populations; subcluster F (brown) contains populations of XLS mountains except for GL. Populations are represented by pie charts, which also represent assignment probabilities between clusters or subclusters. Population codes, sampling information, and number of analyzed individuals are shown in Supplementary Table S1. D: Plot on left shows functional relationship between ΔK and K; K value corresponding to ΔKmax represents optimal grouping; right plot shows membership assignment of East cluster for K=2. E: Left plot shows membership assignment of all populations for K=2; right plot shows maximum value of ΔK corresponding to K=2. F: Left plot shows membership assignment of West cluster for K=4; right plot represents maximum value of ΔK corresponding to K=4. Black lines separate populations. G: Neighbor-joining tree of L. taliangensis populations based on DA distances. Branch color represents subclusters (A–F) results of STRUCTURE, population code color represents two main clusters.

Mitochondrial DNA (mtDNA) analysis

DNA sequences were edited by MEGA v6.06 (Tamura et al., 2013) and aligned using SeaView v4.7 (Gouy et al., 2010). Each protein-coding DNA sequence could be translated into amino acid sequences and no nucleotide gaps were found. Lasergene v7.1 (DNASTAR, USA) was used to combine gene fragments (cyt b+ND2+COI) for further phylogenetic analysis.

DNASP v5.1 (Librado & Rozas, 2009) was used to generate the haplotype file, and to calculate haplotypic diversity (Hd) and nucleotide diversity (π). Geographic patterns of Hd and π were visualized using the inverse distance weighted (IDW) method in ArcGIS v10.3 (ESRI, USA). Haplotype networks were generated with the TCS network module (Clement et al., 2000) in PopART v1.7 (Leigh & Bryant, 2015). PartitionFinder v2.1.1 (Lanfear et al., 2017) was used to identify the best partition scheme of mitochondrial concatenated genes (cytb+ND2+COI) and to explore the best evolutionary model for each partition based on Bayesian Information Criterion (BIC) and the greedy algorithm. Bayesian inference (BI) analyses were conducted in MrBayes v.3.2.2 (Ronquist et al., 2012) with four Markov chain Monte Carlo (MCMC) chains of 20 000 000 generations, with samples taken every 1 000 generations and the first 25% discarded as burn-in. Maximum-likelihood (ML) analyses were performed in RAxML v7.0.4 (Silvestro & Michalak, 2012) with 1 000 bootstrap replicates for the evaluation of branch confidences. Optimal trees were rooted and visualized using Figtree v1.4.2 (Rambaut & Drummond, 2012).

The two-step method (Hedges & Kumar, 2004) was used to estimate divergence time based on mitochondrial concatenated genes in BEAST v1.8.2 (Drummond & Rambaut, 2007). For the first step, four Salamandridae fossil records were used as primary calibration points (C1–C4, GenBank Nos. in Supplementary Table S4) (Wang et al., 2018) to establish the divergence date of the most recent common ancestor (MRCA) of both L. taliangensis and T. pseudoverrucosus; in the second step, the above divergence date was used for secondary calibration, which yielded the divergence times of the internal clades of L. taliangensis. The dataset resulting from the second step in BEAST was used for ancestral area reconstruction in RASP v3.2 (Yu et al., 2015) with the statistical dispersal-vicariance analysis (S-DIVA) model. Four distributional areas were delimited according to the mountains across the distribution range of L. taliangensis: i.e., DLS, XLS, XXL, and GG. The max areas at each node were set to two.

Neutrality tests (Tajima’s D (Tajima, 1989), Fu’s FS (Fu & Li, 1993)), and mismatch distribution analysis (via sum of squares differences (SSD) and Harpending’s raggedness index (Hrag) (Harpending, 1994; Rogers & Harpending, 1992)) were implemented for all populations and two clades using Arlequin v3.0 (Excoffier & Lischer, 2010). These tests were employed to estimate whether the species has experienced historical population expansion. Under the sudden demographic expansion model, the goodness-of-fit of the observed mismatch distribution was tested via 1 000 bootstrap replicates. The mismatch distribution curve was drawn by DnaSP v5.0 (Librado & Rozas, 2009). Bayesian Skyline Plots (BSPs) were drawn in BEAST v1.8.2 to analyze effective population size fluctuations of all populations and two clades through evolutionary time. Results were displayed using TRACER v1.7 (Rambaut et al., 2018).

Further details on mtDNA analysis are shown in the Supplementary Methods.

Microsatellite analysis

Null alleles and allele dropout of microsatellite data were assessed by MICROCHECKER v2.2.3 (Van Oosterhout et al., 2004). Significant deviations from both the Hardy-Weinberg equilibrium (HWE) and linkage disequilibrium between all pairs of microsatellite loci were identified by Genepop v4.2 (Rousset, 2008). Significance values were Bonferroni corrected (Rice, 1989).

Total number of alleles (TA), number of effective alleles (Ne), observed heterozygosity (Ho), unbiased expected heterozygosity (uHe), and inbreeding coefficients (F) for each population were calculated by GenAlEx v6.5 (Peakall & Smouse, 2012). Allelic richness (Ar) and private allelic richness (pAr) were calculated in ADZE v1.0 (Szpiech et al., 2008) and were corrected for the effects of different sample sizes of each population by the rarefaction approach. Geographic patterns of genetic diversity (uHe and Ar) were visualized by the inverse distance weighted (IDW) method in ArcGIS v10.3.

Population structures were investigated using STRUCTURE v2.3.1 (Pritchard et al., 2000) with Bayesian clustering based on microsatellite loci with a burn-in of 100 000 iterations, followed by 1 000 000 MCMC iterations. Number of genetic clusters (K) was tested from 1 to 8 with 10 replicates for each K value. The criterion for optimal K was the highest value of delta K (ΔK) (Evanno et al., 2005), which was calculated in STRUCTURE HARVEST (Earl & vonHoldt, 2012). The outputs of optimal K with 10 replicates were averaged and calculated using CLUMPP v1.1 (Jakobsson & Rosenberg, 2007). The diagram of the best genetic clusters was drawn by DISTRUCT v1.1 (Rosenberg, 2004). According to the population structure, hierarchical analysis of molecular variance (AMOVA) was performed both in whole populations and two main genetic clusters using Arlequin v3.5.1 with 1 000 permutations. Principal coordinate analysis (PCoA) was conducted by GenAlEx v6.5, and the results were visualized in R v3.6.1 using the software package “ggpolt2”. Furthermore, Nei’s genetic distance (DA; Nei et al., 1983) of the microsatellite data was also calculated and a neighbor-joining (NJ) tree was constructed in POPTREE2 (Takezaki & Tamura, 2014) with 1 000 bootstraps.

Heterozygosity excess, which presents evidence of a population bottleneck, was estimated by the sign-test and Wilcoxon sign-rank test using the stepwise mutation model (SMM) and two-phase model (TPM) (proportion of SMM in TPM=95%, and variance of TPM=12). BOTTLENECK v1.2.02 was used with 1 000 replicates (Piry et al., 1999). To estimate the direction of migration and gene flow between major genetic clusters, mutation-scaled migration rates (M) and effective population sizes (θ) of each cluster were calculated using the Bayesian coalescent approach in MIGRATE-N v3.7.2 (Beerli & Palczewski, 2010). Differences in the inbreeding coefficient (F), mean corrected assignment index (mAI), and variance of corrected assignment index (vAI) were compared between sexes in the whole population and two genetic clusters (Goudet et al., 2002). For this, Fstat v2.9.4 (Goudet, 1995) was used with 1 000 permutations to estimate sex-biased dispersal. F andvAI of the dispersing sex are higher than the philopatric sex, while mAI is lower for the dispersing sex (Goudet et al., 2002).

Additional details on microsatellite analysis are shown in the Supplementary Methods.

Ecological niche modeling

To infer the potential distribution range of L. taliangensis under historical and current climatic conditions and to estimate whether populations have experienced contraction during periods of climatic fluctuations, ENM was applied (Elith et al., 2011). This approach combines current distribution GPS information of the species and environmental variable modeling to project historical distributions.

Current distribution information was obtained from field investigations, as well as records from the Global Biodiversity Information Facility (https://www.gbif.org). Only one GPS site within each pixel (30 arc-seconds, 1 km×1 km) was retained. Moreover, 19 predictor climate variables in four periods were downloaded from WorldClim (https://www.worldclim.org). According to the distribution records, the current distribution region and surrounding 200 m buffer area (N25–31°, E100–105°) were selected as the range of interest. ENM was conducted in MaxEnt v3.4.1 (Phillips et al., 2017) using five bioclimatic variables and maximum entropy method based on optimal parameters. The area under the receiver operating characteristic curve (AUC) (Swets, 1988) was used as the evaluation index of the overall fit of the model (Araújo & Guisan, 2006). The MaxEnt output was converted into the binary classification of either suitable (i.e., 1) or unsuitable (i.e., 0) based on the maximum training sensitivity plus specificity (MTSS) threshold (Liu et al., 2016). To identify the climatic stability area, the distribution prediction rasters of four periods were overlapped by ArcGIS and a raster with classification values ranging from 0 (i.e., not present in any of the four time periods) to 4 (i.e., present during all four periods) was generated.

For more detailed information, please refer to the Supplementary Methods.

Landscape genetics

To assess whether the population structure and patterns of gene flow are related to geographical barriers or environmental differences, IBR (which uses resistance distance instead of Euclidean geographical distance) and IBE were examined. These are often used to study the effects of landscape characteristics on genetic structuring (Atlas & Fu, 2019; Myers et al., 2019; Vasconcellos et al., 2019; Waraniak et al., 2019).

To alleviate the greater subjectivity associated with expert experiences (Myers et al., 2019; Spear et al., 2010), the resistance distance for the IBR model was calculated in CIRCUITSCAPE v4.0 (McRae, 2006; McRae et al., 2016) based on the resistance surfaces derived from the ENM result rasters. The environmental dissimilarity matrix, used for IBE analysis, was calculated from 20 environmental parameters. The pairwise genetic distance matrix (F) of the microsatellite data was calculated in Arlequin v3.5.1. All above distance matrices were standardized by R v3.6.1 using the “scale” function. The variance inflation factor (VIF) was calculated before regression analysis to avoid multicollinearity of independent variables. The multiple matrix regression with randomization (MMRR) function was employed to evaluate the significance of effects of multiple independent factors (environmental dissimilarity and resistance distance in this study) on genetic divergence (Wang, 2013), conducted in R v3.6.1 with 10 000 permutations. The backward-stepwise method was used to select the optimal model. In addition, significance of the rejected variables was also tested.

Further details on resistance surface construction and environmental dissimilarity matrix calculation are presented in the Supplementary Methods.

Dispersal route analysis

To evaluate the influence of landscape features on the connectivity of dispersal routes of L. taliangensis, least-cost corridors (LCCs) and least-cost paths (LCPs) were calculated using SDMTOOLBOX v1.1c (Brown, 2014) in ArcGIS v10.3. The resistance raster used in IBR analysis was used as an analysis raster and sampling site coordinates were entered in decimal degrees. The LCPs were weighted by resistance values and classified into three categories using the “percentage of LCP” method (cutoff values: 1%, 2%, and 5%); finally, an LCC dispersal network was created.

RESULTS

Mitochondrial phylogeography

After manual sequence correction, alignment, and trimming of ends, concatenated gene sequences with a total length of 3 255 bp were obtained with three gene fragments (i.e., ND2 (945 bp), cyt b (762 bp), and COΙ (1 548 bp)) from 407 samples. No gaps were found, and all sequences could be correctly translated into amino acid sequences based on vertebrate genetic codes. The sequence data of three mtDNA haplotypes were deposited in GenBank under accession Nos.: cyt b (MZ368891–MZ368957), ND2 (MZ368958–MZ369024), and COI (MZ369025–MZ369091).

For the combined gene sequence, 156 mutation sites were found, and 67 haplotypes were defined (i.e., H1–H67), including 35 private haplotypes belonging to six populations (Supplementary Table S5). The haplotype network showed that adjacent clades or subclades were separated by at least seven mutational steps. The non-single star-shaped haplotype network suggested that the population has experienced multiple local expansions (Supplementary Figure S1). Populations Yuexi (YX) and Chuhongjue (CHJ) exhibited the highest population haplotype diversity and highest nucleotide diversity, respectively (Supplementary Table S6). The North clade had a much higher Hd and π than the South clade (Figure 2; Supplementary Table S7).

Figure 2

Spatial interpolation of genetic diversity parameters among populations of Liangshantriton taliangensis

A: Allelic richness (Ar); B: Unbiased expected heterozygosity (uHe); C: Haplotype diversity (Hd); D: Nucleotide diversity (π). Warmer colors indicate higher genetic diversity. Black dots represent sampling sites. Further information on population codes is given in Supplementary Table S1.

The best substitution models (i.e., K80+I for the 1st COI and cyt b codon partition, HKY+I for the 2nd COI, cyt b and ND2 codon partition, GTR+G for the 3rd COI, cyt b and ND2 codon partition, HKY+G for 1st ND2 codon partition) were identified using PartitionFinder v2.1.1. Both BI and ML phylogenetic trees clearly divided L. taliangensis into two monophyletic lineages (i.e., North and South clades, BI/ML: 100%/100) (Supplementary Figure S1). Both clades had mostly allopatric distributions, with the South clade concentrated in the DLS, and the North clade concentrated in the XLS and XXL-GG (Supplementary Figure S1). The North clade could also be divided into seven lineages (subclades B–H) , although the relationships between specific subclades could not be completely resolved (Supplementary Figure S1). Subclades B, G, and H were from XXL-GG and Subclades C–E were from XLS. The Xinmin (XM) population from GG was isolated as Subclade H. Subclade F consisted of four populations (i.e., Jinghuahu (JHH), Pusagang (PSG), Zima (ZM), and GL) from XXL and XLS (Supplementary Figure S1 and Table S5).

Divergence time results indicated that L. taliangensis diverged from its sister species T. pseudoverrucosus during the early Pleistocene period, approximately 2.44 million years ago (Ma) (95% highest posterior density (HPD): 1.65–3.28 Ma) (Supplementary Figure S2). The most recent common ancestor (MRCA) of L. taliangensis was estimated at 1.49 Ma (95% HPD: 1.07–1.97 Ma) and the first vicariance event occurred simultaneously (Figure 3). The seven subclades formed before the Guxiang Glaciation (0.3–0.13 Ma) and underwent four dispersal and three vicariance events. The Zhuma (ZUM) and GL populations were the first colony for XXL and transit point between XLS and XXL, respectively. In addition, two dispersal events (both from DLS to XLS) and two vicariance events occurred within the South clade during the Guxiang glacial-interglacial and last glacial-interglacial cycles (Figure 3). The LS was the most likely original distribution area ofL. taliangensis (node 133), with a probability of 66.1% (Figure 3), with the DLS geographically closest to the T. pseudoverrucosus outgroup. The BSP for all populations and the North clade suggested that the species experienced a strong decrease in effective population size after the Last Interglacial (LIG), followed by a very recent increase after the Last Glacial Maximum (LGM) (Figure 4A, C). However, the BSP for the South clade confirmed general stabilization during the LGM and accelerated expansion after the LGM (Figure 4B). Neutrality and mismatch distribution analyses basically supported these results (Supplementary Figure S3 and Table S7).

Figure 3

Ancestor distribution areas of Liangshantriton taliangensis inferred from S-DIVA analyses based on mitochondrial data

Figure 4

Bayesian skyline plots (BSP) showing variation in effective population sizes through time

Sectorial annular represents proportions of ancestral range and different ancestral ranges are denoted by different colors. Blue arrows represent direction of dispersal event, blue circles represent dispersal events, and red circles represent vicariance events. Number in sectorial annular is node ID. Divergent times (Ma; based on mitochondrial concatenated genes) are annotated beside main nodes. Numbers in square brackets represent 95% HPD of divergence time. Clades and subclades are named on the right.

A–C: For all populations (A), South clade (B), and North clade (C) of Liangshantriton taliangensis based on mitochondrial concatenated gene sequences. Blue area represents 95% confidence interval of effective population size. Ranges within red dotted lines represent LGM and ranges within blue dotted lines represent LIG.

Genetic structure and dispersal pattern based on microsatellite data

No significant linkage disequilibrium was detected at any loci after Bonferroni correction. Except for homozygous loci in the population, which could not be evaluated for HWE, two loci showed significant departures from HWE in one to six populations following Bonferroni correction. DLY32 indicated significant departure from HWE in one population (i.e., Sanhekou (SHK)), while DLY96 showed significant departure in five populations (i.e., XM, ZUM, Gongyihai (GYH), ZM, and PSG). However, this may be the result of the population genetic structure. In addition, deviations from HWE were not consistent across populations, with most loci in the population conforming to HWE. The microsatellite genetic diversity of each population was estimated (Supplementary Table S8). Ho (mean=0.488) was lower than uHe (mean=0.506), indicating that the population suffered from heterozygosity deficiency, and the mean value was 0.308 for Ar (2.182–4.389). The higher genetic diversity was mostly concentrated in LS, together with the ZUM population in XXL (Figure 2). The positive F values (mean F=0.022>0) suggested that all populations had an inbreeding tendency. The results of population genetic distance (F) showed that there was significant genetic divergence (P<0.001), and the average genetic distance was 0.342 (0.033 (JHH, PSG)–0.598 (ZM, Bingtu (BT))) (Supplementary Table S9).

Structural analysis showed that ΔK reached the highest peak at K=2, indicating that all individuals could be divided into two genetic clusters (Figure 1A, E). The red cluster was defined as the West cluster (XXL-GG) and the blue cluster was defined as the East cluster (LS) (Figure 1A). Relatively high admixture was detected where the genetic clusters were in contact, particularly in the ZUM population (assignment proportion: 68.7% from XXL-GG and 31.3% from LS). Furthermore, four subclusters were revealed in the West cluster (Figure 1B, F) and K=2 was the optimal grouping for the East cluster (Figure 1C, D). Genetic diversity of the East cluster was higher than that of the West cluster (Supplementary Table S8). Similarly, PCoA suggested that all populations could be divided into two genetic clusters (i.e., East and West cluster), which did not overlap on the PCoA1 axis (Supplementary Figure S4). The NJ tree (inferred from DA distances) also identified two clusters and six subclusters (Figure 1G). Subclusters C and F showed the highest genetic divergence (Figure 1G). AMOVA showed that genetic variation was slightly lower between two clusters than among populations within clusters (Supplementary Table S10). Within the West cluster, genetic variation was much higher among groups than among populations within groups. However, within the East cluster, genetic variation was lower among groups than among populations within groups and no significant inbreeding was found (Supplementary Table S10). The TPM model, which is suitable for microsatellite data (Di Rienzo et al., 1994), indicated that L. taliangensis was not subjected to a bottleneck effect (Supplementary Table S11).

The migration patterns among genetic clusters were from XXL-GG to LS, with the highest Bezier approximation score of −71 223.03 lmL (Figure 1; Supplementary Table S12). Under this model, the mutation-scaled population size of XXL-GG (ΘXXL-GG=5.97, 95% CI: 2.27–9.47) was almost two times higher than that of LS (ΘLS=2.99, 95% CI: 0.00–6.40) and had a very low M value (i.e., immigration rate/mutation rate) (MX-L=0.467<1) for one-way gene flow. These results indicated thatL. taliangensis experienced very limited unidirectional gene flow from XXL-GG to LS.

F,mAI, and vAI were compared between females and males over the whole population. Results showed that F was significantly higher in males than in females, indicating that males were the dispersing sex, while vAI was significantly higher for females than for males, indicating that females were the dispersing sex. mAI showed no significant differences between the sexes. Thus, dispersal sex was not fixed within the whole population (Table 1). We checked the two main clusters, respectively. Only mAI showed a significant difference, with F and vAI showing non-significant sex differences in the two clusters. For the West cluster, the mAI values were significantly lower in males than in females, suggesting male-biased dispersal; for the East cluster, the mAI values were significantly lower in females than in males, indicating that females were the dominant dispersal sex (Table 1).

Table 1

Differences in inbreeding coefficient (F), mean of corrected assignment index (mAI), and variance of corrected assignment index (vAI) between males and females of L. taliangensis from total range, West cluster, and East cluster

	N	F IS	mAI C	vAI C
N: Sample size. *: 0.01<P<0.05.
ALL
Females	83	–0.023	0.536	82.169
Males	324	0.043	–0.137	50.659
P		0.014*	0.533	0.028*
West cluster (XXL-GG)
Females	49	–0.014	1.645	38.690
Males	152	0.072	–0.530	57.506
P		0.061	0.046*	0.120
East cluster (LS)
Females	34	–0.023	–3.949	86.440
Males	172	0.024	0.781	54.285
P		0.247	0.017*	0.133

Changes in distribution area during glaciations

Linear and quadratic features and a regularization multiplier of 0.5 (LQ 0.5) were optimal model parameters, and the very high AUC values (0.963±0.018) indicated that this model had high predictive power. The MTSS threshold was 0.2064. The ENM predicted wider areas of suitable L. taliangensis habitat in the tropical and subtropical regions of the southeastern HMR during the LIG than during the LGM, mid-Holocene (MH), or the present (Figure 5). Based on the distribution predictions from the three general circulation models (GCMs), the potentially suitable distribution areas ofL. taliangensis shrank during the dry and cold LGM (Figure 5B; Supplementary Figure S5). The predicted distribution range during the MH was similar to the current range but much wider (Figure 5C, D; Supplementary Figure S6). Overlap of the predictions in four periods showed that the climatic stability zone was mainly distributed in the current distribution range (Figure 5E). Isothermality (Bio03), minimum temperature of coldest month (Bio06), annual precipitation (Bio12), and precipitation of driest month (Bio14) were high in contribution and permutation importance. The cumulative contribution rate was 97.59% (Supplementary Table S13). The most critical climatic predictor in the model was minimum temperature of coldest month (Bio06), with a contribution rate of 36.93% (Supplementary Table S13).

Figure 5

Species distribution range of Liangshantriton taliangensis under four different periods and overlapping map yielding historical climatic stability areas

A: Last Interglacial (LIG). B: Last Glacial Maximum (LGM). Projection of LGM was the average of three general circulation models (GCMs) (CCSM4, MIROC-ESM, and MPI-ESM-P) (see Supplementary Figure S5). C: mid-Holocene (MH). Projection of MH was the average of four GCMs (CCSM4, MIROC-ESM, MPI-ESM-P, and BCC-CSM1-1) (see Supplementary Figure S6). D: Present. Based on MTSS, which represents maximum training sensitivity plus specificity threshold, projection distribution area was divided into areas with low, medium, and high suitability, denoted in green, yellow, and red, respectively. E: Showing historical climatic stability areas, from presence in only one period to presence during all four periods. Two major genetic clusters are represented by hollow triangles and circles.

Critical landscape factors affecting population structure

The ENM used for IBR analysis showed good performance in predicting the potential distribution of L. taliangensis (AUC=0.936) and MTSS was 0.2301 (Supplementary Figure S7). The VIF did not show significant multicollinearity between resistance distance and environmental dissimilarity (√VIF<2). A multiple matrix regression model combining resistance distance and environmental dissimilarity was generated via MMRR. The model indicated that genetic distance (F) of all populations (Supplementary Table S9) was significantly positively correlated with resistance distance (R2=0.115, P=0.013) and all environmental variables were excluded from the final model (Figure 6A; Supplementary Table S14). However, the removed factors were used for repeat MMRR analysis, which showed that only PC1 (temperature-related variables, Supplementary Table S15) was significantly, though weakly, positively correlated with genetic distance (R2=0.060, P=0.019) (Figure 6B; Supplementary Table S14). For the West cluster, only resistance distance (IBR) strongly explained genetic divergence (R2=0.938, P=0.001) (Supplementary Figure S8 and Table S14). No landscape characteristics were associated with genetic divergence in the East cluster (Supplementary Table S14). Evaluations of variable contribution and permutation importance identified Bio14 (precipitation of driest month) and altitude as the most important factors affecting species distribution and dispersal, while normalized difference vegetation index (NDVI) and Bio12 (annual precipitation) were the next most significant variables (Supplementary Table S16).

Figure 6

Multiple matrix regression with randomization (MMRR) plots for isolation-by-resistance (IBR) and isolation-by-environment (IBE) as well as construction of least cost corridors (LCCs) forLiangshantriton taliangensis based on resistance raster

A: Correlations between genetic distance (F) and resistance distance. B: Correlations between genetic distance (F) and PC1 environmental distance. C: LCCs among 16 sampled populations of L. taliangensis. Warmer colors indicate higher connectivity and lower dispersal cost. Sampling sites are represented by black dots.

Dispersal route analyses

According to the LCC map, neighboring populations had better connectivity than distantly separated populations (Figure 6C). There were few corridors with low resistance between DLS and both XLS and XXL-GG (Figure 6C). The dispersal corridor between the XM population and other populations showed high geographic and environmental resistance, implying that individuals from the XM population and other populations were not easy dispersal each others. Of note, a narrow dispersal corridor with low resistance through the GL and Niuri River (tributary of the Dadu River) between XXL and XLS indicated the possibility of gene flow between populations of both mountains.

DISCUSSION

Regional uplift and Early-Middle Pleistocene transition (EMPT) drove L. taliangensis formation and divergence

The causes of species evolution and diversification are diverse and complex, and generally include both biotic and abiotic factors or a combination of both, e.g., climate change, geographic barriers, and inter- and intra-specific interactions (Bouchenak-Khelladi et al., 2015; Doebeli & Dieckmann, 2003). In this study, we discovered that L. taliangensis originated from the LS mountains, separated from its sister species (T. pseudoverrucosus) during the Pleistocene (~2.44 Ma) (Supplementary Figure S2), and dispersed to higher latitudes. This is mostly consistent with the divergence times and radiation direction of Tylototriton s.l. (Wang et al., 2018), although the divergence time estimated by secondary calibration should be interpreted with caution (Hedges & Kumar, 2004; McCormack et al., 2011).

The QTP uprising promoted the formation of the East Asian monsoon, creating a warm, moist, and vegetation-rich habitat on the southeastern edge of the HMR, which was suitable for newts (Wang et al., 2018). Climatic shifts and increased niche vacancies drove the radiation of ancestral populations of L. taliangensis from low latitudes to high latitudes. The South and North clades separated ~1.49 Ma (Figure 3), with the uplift of the southeastern HMR. The first vicariance event was detected during this period. Although details on the timing and pattern of uplift remain controversial, consensus has been achieved that the uplift of the HMR happened recently and rapidly (Xing & Ree, 2017). Speciation and divergence of amphibians are presumed to be the result of ecological isolation caused by uneven uplift of the QTP (Hofmann et al., 2017; Wu et al., 2020). About 2.5–1.2 Ma, the Xigeda Formation strata on the southeastern margin of the QTP (widely distributed between the Dadu River, Anning River, and Jinsha River valley in southwestern Sichuan) gradually began to uplift and complete the accretionary orogenic process (Jiang & Wu, 1998). Consequently, the formation of species and divergence of L. taliangensis clades may be related to the East Asian monsoon climate and regional uplift events (e.g., Xigeda Formation strata). The LS mountains (DLS and XLS) showed the highest probability as the ancestral area for L. taliangensis. The existence of multiple ancestral areas suggests that the divergence patterns of the main clades may be consistent with the fragmentation model, where isolation and fragmentation occurred over a wide range of ancestral distributions.

The North clade of L. taliangensis underwent phylogenetic divergence around the EMPT, although the phylogenetic relationships among subclades remain unclear. The EMPT (~1.4–0.4 Ma, Head & Gibbard, 2015) was characterized by a gradual increase in the amplitude of climatic oscillations and in global ice volume, with a shift from 41 thousand year (ka) glacial cycles to 100 ka cycles since ~1.2 Ma (Clark et al., 2006; Head & Gibbard, 2015). From ~1.25 Ma onwards, the East Asian summer monsoon during the interglacial periods intensified (Sun et al., 2010), which further promoted the dispersal of L. taliangensis from XLS to XXL-GG. At ~0.9 Ma, a dispersal and vicariance event of L. taliangensis was recorded, simultaneously. This period, known as the “0.9 Ma event”, was the beginning of the first large-scale expansion of continental ice sheets and was interrupted by the weak interglacial period (Marine Isotope Stage (MIS) 23 ) (Clark et al., 2006; Head & Gibbard, 2015). A total of three dispersal and two vicariance events occurred during the early phase of the EMPT, with divergence of the XLS population from the XXL-GG population. The Naynayxungla Glaciation (0.78–0.5 Ma) had a profound influence on the formation of the internal lineage of the XXL-GG, and one dispersal and one vicariance event occurred during this period. The XM populations diverged into a monophyletic lineage at ~0.77 Ma. This time point, called the Matuyama-Brunhes boundary (773 ka), is the primary guide for the Lower-Middle Pleistocene Subseries boundary (Head & Gibbard, 2015). During this period, the global temperature decreased significantly, ice sheets expanded, and glacial cycles were dominated by 100 ka periodicity (Clark et al., 1999, 2006; Head & Gibbard, 2015). Furthermore, the QTP had a glacier coverage of close to 500 000 km2 (Zheng et al., 2002). The Shaluli Mountain region and Daocheng Kuzhaori, which are relatively close geographically to the current distribution area of L. taliangensis, were the main areas and more strongly influenced by QTP glaciation (Wang et al., 2006; Xu et al., 2005). A dispersal event from XXL to GL occurred ~0.49 Ma during the great interglacial period, which was related to the strong summer monsoon throughout the Northern Hemisphere (Figure 3). The overall dispersal direction was from low to high latitudes, and ZUM was identified as the first colony of XXL (Figure 3). Subclade divergence mainly supported a stepping-stone colonization mode, with dispersal from one mountain to another, followed by vicariance. GL was an important stepping-stone between XLS and XXL (Figure 3).

Frequent vicariance and dispersal of populations due to Pleistocene climatic fluctuations may be a common cause of lineage divergence in many rapidly diverging species (Wang et al., 2017, 2018; Ye et al., 2019), including L. taliangensis. The formation and differentiation of L. taliangensis also provide further biological reference for the timing controversy associated with orogenesis.

Multiple glacial refugia patterns respond to LGM

Refugia tend to maintain moderate temperature and precipitation levels, and are thus commonly identified as climatically stable regions, which enabled species to survive during the LGM (Tang et al., 2018; Tzedakis et al., 2002). Different from Europe and North America, there was no unified ice sheet covering the QTP in China and many ice-free low-altitude areas could be found on its southeastern margin (Owen et al., 2008). The southeastern HMR was a large glacial refugium for endemic species and its complex topographic features allowed for lineage maintenance in “refugia within refugia” (Xie et al., 2019; Yu et al., 2013). In this study, the current distribution area of L. taliangensis represents such a climatically stable zone, mainly located in the Dadu River valley and the LS and XXL intersection area (Figure 5E). Liangshantriton taliangensis showed multiple glacial refugia patterns in response to the LGM, which is common in species with weak dispersal capacity in southwest China (Liu et al., 2015; Wang et al., 2017). The single refugium hypothesis suggests that genetic diversity decreases from glacial refugia to surrounding areas, and refugia usually harbor private haplotypes and alleles (Hewitt, 2000). The following glacial refugia forL. taliangensis were identified (Figures 2, 5E; Supplementary Tables S6, S8): one in GG (XM), two potential major glacial refugia in XXL (ZUM and JHH), two in DLS (YX and Qiliba (QLB)), and two main glacial refugia in XLS (GL and CHJ). The response patterns of population dynamics to Pleistocene climatic fluctuations vary among species and distribution ranges (Lu et al., 2012; Qu et al., 2010). Here, two major lineages (i.e., South and North clades) displayed diverse population dynamics (Figure 4). The South clade population was located at relatively low latitudes and on medium-high altitude mountains below the snow line. Therefore, L. taliangensis in this region was probably rarely affected by glaciation. The GG mountains are considered as the glaciation center of the HMR, where mountains above the snow line became a source of cold, the influence of which extends from north to south (Li et al., 1991). The North clade (i.e., GG, XXL, and XLS) population was distributed at higher latitudes and adjacent to the GG mountains. Thus, L. taliangensis in this range was subjected to glacial influences.

Mitonuclear discordance

Amphibians, especially frogs, commonly show mitonuclear discordance; however, it is rare in salamanders (Chan & Levin, 2005; Bryson et al., 2014; Toews & Brelsford, 2012). In general, mitonuclear discordance results from female-biased dispersal, incomplete lineage sorting, hybridization/introgression in cases of secondary contact, and selective sweep of mtDNA haplotypes (Garrick et al., 2019; Lucati et al., 2020; Toews & Brelsford, 2012). In this study, female-biased dispersal was found in the East cluster (Table 1), and introgression by secondary contact was identified in several regions and at different times (Figure 3). Thus, we speculated that these two factors may be the main reasons for the observed mitonuclear discordance.

Inferred from mtDNA, L. taliangensis could be divided into two well-supported genetic lineages, i.e., the South and North clades (separated at ~N28.6°). In contrast, two main genetic clusters of L. taliangensis were recognized by the microsatellite data, i.e., West and East clusters, which exhibited a more explicit biogeographic pattern and were separated by the XXL mountains (Figure 1), similar to the sympatric red pandas (Ailurus fulgens F. Cuvier, 1825) (Hu et al., 2011). Although contrasting divergence patterns, associated with the main genetic structures in the mitochondrial and nuclear genotypic data, were detected in the current study, fine-scale structures within these main genetic clusters showed strong concordant patterns.

Within the East cluster, two subclusters (i.e., E and F) were recognized. Subcluster E mainly contained populations from DLS (YX, QLB, and Wuke (WK)), and partly corresponded to the mitochondrial South clade at the base. Ancestral range reconstruction analysis identified this as a plausible independent biogeographic area. Therefore, if populations from DLS were to be accepted as an independent genetic cluster, the BT population and a few individuals from the CHJ population, with mitochondrial haplotypes nested within the South clade, may be the result of introgression due to secondary contact between the South and North clades. Two secondary contacts may have occurred between the South and North clades based on ancestral range reconstruction analysis. Combined with haplotype information (Supplementary Table S5), the first introgression may have occurred ~0.34 Ma, with dispersal from the QLB population (DLS) to the BT and CHJ populations (XLS); the second introgression was a dispersal from the YX population (DLS) to the CHJ population (XLS), which occurred during the LIG (Figure 3). Sex-biased dispersal suggested that the LS population was dominated by female dispersal, indicating that females drove asymmetrical mtDNA introgression (Table 1). The CHJ population may represent an important secondary contact zone between DLS and XLS.

It is worth noting that neither mitochondrial nor microsatellite data showed significant gene flow between GL population and other populations in the XLS (Figure 1; Supplementary Figure S1). Although dispersal route analysis revealed a narrow low-resistance dispersal corridor between GL and both SLP and CHJ populations, and weak nuclear gene mixing between them was detected by microsatellites, the limited barrier effect of Ma'an Mountains (altitude: 4288 m a.s.l.) cannot be ignored. Mitochondrial data suggested secondary contact between the GL and XXL populations at ~0.49 Ma and isolation from the XXL population during the Guxiang Glaciation (0.3–0.13 Ma), with mitochondria of the GL population partially replaced by that of the XXL population (Figure 3). The good landscape connectivity between the GL and XXL populations (especially GYH and JHH) and migration analysis both suggested the possibility of limited gene flow from XXL to XLS (Figures 1, 6C). GL may represent an important secondary contact zone between XXL and XLS. The GL population belonged to the same nuclear genotypic cluster as populations from DLS. However, this pattern contrasted with the mitochondrial phylogeographic patterns and dispersal route analysis, which indicated high dispersal resistance between the DLS and GL populations and a low potential for gene flow (Figure 6C). These findings suggest that the distinct genetic differences of the GL population from other populations may be diminished by size homoplasy of microsatellites or strong selection pressure on sequence-based loci, which often occur in populations with large effective population size (Epperson, 2005). However, further studies are needed to confirm this speculation using more robust genetic and morphological data.

All populations from XLS (except for GL) belonged to subcluster F, which was mainly associated with the two mitochondrial subclades (C and D). This divergence in the mitochondrial clades occurred during the Pleistocene, as tracers of Quaternary glacial events, indicating that repeated expansions and contractions of the population ranges was accompanied by lineage sorting. However, the homogenization of nuclear genotypic data and mitochondrial introgressions between subclades C and D may reflect recent gene flow, resulting from population expansion after the Quaternary. Within the XXL-GG, two subclusters corresponded to two mitochondrial subclades. Apparent mitonuclear concordances were reflected by nuclear subcluster A and mitochondrial subclade H (XM population) as well as nuclear subcluster D and mitochondrial subclade B (from the ZUM population). Discordance in the distribution of nuclear genotypes with mitochondrial clades mainly manifested in the BT, CHJ, GL, and GYH populations. The GYH population showed a distinct nuclear genotype but was characterized by mitochondrial subclade G haplotypes. The geographical approximation between GYH and other populations of subclade G and the high landscape connectivity in this region suggested frequent gene flow (Figure 6C); thus, mitochondrial capture from neighboring populations may explain the above (Hausdorf et al., 2011).

Dramatic climatic fluctuations and scattered glacial refugia can promote the occurrence of frequent population expansion and isolation, which may be an important cause of mitonuclear discordance. These effects provide opportunities for gene admixture and fusion, and traces of asymmetric gene flow at distribution boundaries are identifiable for longer in mitochondrial than in nuclear genes (Currat et al., 2008; Dufresnes et al., 2020; Excoffier et al., 2009). Here, GL and CHJ in XLS were important glacial refugia and major secondary contact zones for L. taliangensis, they also distributed at the geographic edge of the intersection of two different genetic clades/clusters.

Resistance distance plays a dominant role in population structure

Resistance distance was a main reason for the population genetic structure of L. taliangensis (Figure 6A; Supplementary Table S14), suggesting that strong gene flow may have occurred in adjacent regions and regions with high geographic connectivity. The barrier effect of geographic distance is ubiquitous among organisms, especially in amphibians due to their low dispersal ability (Bani et al., 2015; He et al., 2021; Waraniak et al., 2019). In addition, the effect of temperature difference was weak but significant when environmental dissimilarity was analyzed alone (Figure 6B; Supplementary Table S14). The West cluster, with its finer population structure and large genetic variation, was strongly influenced by geographic resistance (Supplementary Figure S8 and Table S14), while the genetic structure of the East cluster was not significantly influenced by landscape features (Supplementary Table S14).

Precipitation of the driest month and altitude were primary factors influencing gene flow between populations, while NDVI and annual precipitation were secondary factors (Supplementary Table S16). Precipitation is crucial for the survival and dispersal of species (especially during the dry season), not only for amphibians but also for reptiles, even though they are perfectly adapted to terrestrial life (Myers et al., 2019). Heavy precipitation contributes to the establishment of temporary connections between adjacent water systems, which facilitate salamander dispersal (Wu et al., 2013). A 70% reduction in forest connectivity can diminish dispersal and further reduce gene flow among amphibians (Campos et al., 2020; Curtis & Taylor 2004; Gibbs, 1998). Salamanders have a particularly high demand for vegetation, not only as terrestrial shelters (Belasen et al., 2013; Otto et al., 2014) but also as microenvironments created by the interaction between precipitation and vegetation (Escoriza & Hassine, 2014). The barrier effect on gene flow, imposed by elevational differences, has been demonstrated in many amphibians, and is particularly important for high-altitude and pond-breeding species (Atlas & Fu, 2019; He et al., 2021; Sánchez-Montes et al., 2018). The north-south running XXL mountains, with a highest peak of 4 791 m a.s.l., may have driven divergence among the East and West genetic clusters ofL. taliangensis. However, these mountains were not absolute barriers (Sánchez-Montes et al., 2018) and limited gene flow could still be detected between the two genetic clusters (Figures 1, 6C). The elevational range of XXL-GG mountains is about 780–5 793 m a.s.l., consisting of steep mountains that are high in the northwest and low in the southeast. Such topographic features may have fragmented the distribution of the West cluster populations. However, the LS mountains consist of the northern XLS with steep topography and altitude range of ~500–4 500 m a.s.l. as well as the southern DLS with a gentle topography and altitude range of about 2 000–3 500 m a.s.l. (Yao et al., 2017). No landscape features that critically affected the genetic structure of the East cluster populations were found. Whether the relationship between population genetic differences and landscape features in the East cluster was masked by variations in the resistance effects of different topographies on gene flow, the size homoplasy of microsatellites, and strong selection pressure on sequence-based loci still need to be further clarified. The population differentiation ofL. taliangensis showed weak correlation with temperature (IBE), especially during the dry and cold seasons, which probably reflects differential adaptation and natural selection to the local climate (Foll & Gaggiotti, 2006; Zhang et al., 2016).

Conservation implications

Research on genetic diversity and population structure provides an important reference and foundation for the development of conservation measures (Ouborg, 2010).

Based on mitochondrial and microsatellite data, the overall genetic diversity of L. taliangensis was high and did not show a bottleneck event (Supplementary Tables S6, S8, S11). None of the sampled L. taliangensis populations currently suffer from genetic deterioration. Although evolutionarily significant units (ESUs) could not be clearly defined (according to Moritz’s strict criteria for ESUs, which require significant differentiation in both mitochondrial and nuclear genes (Moritz, 1994)), two major genetic clusters (i.e., West and East clusters) could be separated as two management units. As critical regions for the persistence and evolution of biodiversity, glacial refugia have conservation priority. Former glacial refugia, including XM, ZUM, JHH, GL, CHJ, YX, and QLB, harbor rich genetic diversity. GL and CHJ are crucial secondary contact zones, and their protection should be prioritized. In addition, attention should be paid to maintaining interpopulation connectivity, especially between XLS and XXL (Figure 6C). Although precipitation, altitude, and vegetation affect the construction of natural corridors (Supplementary Table S16), human population expansion, agricultural development, deforestation, and related problems may put these corridors at risk of disappearing, especially in XXL (Hu et al., 2011). Therefore, field surveys should be conducted for the establishment of alternative corridors while maintaining those that exist.

CONCLUSIONS

Based on multiple genetic markers (i.e., mitochondrial concatenated genes (cyt b+ND2+COI) and 11 microsatellite loci), the genetic divergence and driving factors ofL. taliangensis, a vulnerable amphibian species endemic in the southeastern HMR, were analyzed. Regional (Xigeda Formation strata) uplift events, the East Asian monsoon, glacier development during the EMPT, multiple in situ refugia patterns during glacial-interglacial cycles, and geographic resistance (precipitation, altitude, and vegetational cover) were identified as the main driving forces for the rapid divergence and population genetic structure formation of this species. Furthermore, secondary contact and/or female-biased dispersal were found to be the main reasons for the observed mitonuclear discordance. The two main microsatellite markers confirmed genetic clusters should be considered as independent management units. Seven major glacial refugia as well as crucial natural corridors between XXL and XLS should also be protected as a priority. This study improves our understanding of the contribution of historical and contemporary environmental factors to the genetic divergence of species in the HMR and provides valuable knowledge for the conservation and genetic management of L. taliangensis.

Supplementary data to this article can be found online.

Click here for additional data file.

COMPETING INTERESTS

The authors declare that they have no competing interests.

AUTHORS’ CONTRIBUTIONS

F.X., J.P.J., C.L., and Z.B.S. designed the study. X.X.S, Y.M.H., M.Y.C., and G.C.S. collected and processed samples. X.X.S., Y.M.H., G.C.S., B.W., and X.Q.L. analyzed the data. X.X.S. and Y.M.H. interpreted the data and wrote the manuscript. F.X. and J.P.J. improved the manuscript. All authors read and approved the final version of the manuscript.