Skip to main content
Download PDF
  • Research article
  • Open access
  • Published:

Dynamics behind disjunct distribution, hotspot-edge refugia, and discordant RADseq/mtDNA variability: insights from the Emei mustache toad

BMC Evolutionary Biology volume 20, Article number: 111 (2020) Cite this article

Abstract

Background

The distribution of genetic diversity and the underlying processes are important for conservation planning but are unknown for most species and have not been well studied in many regions. In East Asia, the Sichuan Basin and surrounding mountains constitute an understudied region that exhibits a “ring” of high species richness overlapping the eastern edge of the global biodiversity hotspot Mountains of Southwest China. We examine the distributional history and genetic diversification of the Emei mustache toad Leptobrachium boringii, a typical “ring” element characterized by disjunct ranges in the mountains, by integrating time-calibrated gene tree, genetic variability, individual-level clustering, inference of population splitting and mixing from allele frequencies, and paleoclimatic suitability modeling.

Results

The results reveal extensive range dynamics, including secondary contact after long-term isolation via westward dispersal accompanied by variability loss. They allow the proposal of a model that combines recurrent contractions caused by Quaternary climatic changes and some failed expansions under suitable conditions for explaining the shared disjunct distribution pattern. Providing exceptional low-elevation habitats in the hotspot area, the eastern edge harbors both long-term refugial and young immigrant populations. This finding and a synthesis of evidence from other taxa demonstrate that a certain contributor to biodiversity, one that preserves and receives low-elevation elements of the east in this case, can be significant for only a particular part of a hotspot. By clarifying the low variability of these refugial populations, we show that discordant mitochondrial estimates of diversity can be obtained for populations that experienced admixture, which would have unlikely left proportional immigrant alleles for each locus.

Conclusions

Dispersal after long-term isolation can explain much of the spatial distribution of genetic diversity in this species, while secondary contact and long-term persistence do not guarantee a large variation. The model for the formation of disjunct ranges may apply to many other taxa isolated in the mountains surrounding the Sichuan Basin. Furthermore, this study provides insights into the heterogeneous nature of hotspots and discordant variability obtained from genome-wide and mitochondrial data.

Background

Understanding the spatial distribution of genetic diversity within species and the underlying processes is important for conservation [1,2,3]. Quaternary climatic oscillations have deeply affected the distribution and genetic diversity of a variety of organisms, especially in the Northern Hemisphere [4, 5]. The resulting relationship between distributional history and diversity is often conditional. For instance, long-term refugia may maintain populations with large effective sizes and high genetic variation, yet drift and inbreeding can cause a reduction in diversity when the refugial areas are small [4, 6]. Alternatively, range expansions are usually characterized by reduced diversity due to the founder effect and the fact that only a portion of all diversity is involved in the expansion [7, 8], while admixture results in increased levels of variation in areas of secondary contact [9]. Such complexities may complicate the prediction of genetic diversity based on the duration of occupation or other measures of unevenness across the geographic range of a species.

Regional biogeographic studies have often identified different major processes that shape the patterns of intraspecific genetic variability, for example, evolutionary melting pots [9, 10], historical climate stability [11,12,13], and reduction of diversity during dispersal [14, 15]. In relatively well-studied areas, commonalities may emerge regarding the relationship between the patterns of the distributional history and genetic diversity. On the Iberian Peninsula, for different plants and animals, genetically highly divergent populations have persisted in different glacial refugia throughout the Pleistocene; in some of these species, postglacial expansion has occurred and led to less diverse populations at northern latitudes [16, 17]. Such generalizations not only facilitate the understanding of biodiversity patterns but may also provide the empirical basis for predicting the genetic diversity pattern of a species from the inferred distributional history in corresponding regions [11, 18, 19]. Currently, for most organisms, genetic diversity characterization is missing or incomplete. Numerous phylogeographic studies conducted before the widespread availability of high-throughput sequencing used a limited number of loci (often one), and obtained measures of genetic diversity, such as heterozygosity, that depend strongly on a large number of loci [20, 21]. In understudied, species-rich regions, performing phylogeographic studies with comprehensive sampling should be one of the first steps in understanding and conserving biodiversity.

The Sichuan Basin and surrounding mountains (Fig. 1) in southwestern China constitute one such region. The topography of these basin-mountain systems was mainly shaped by tectonic activity before the Late Pliocene (3.6–2.6 million years ago; mya) [22,23,24,25]. During the Pleistocene, this region remained unglaciated except for a few marginal high-elevation areas in the west [26, 27]. Within this region, many organisms are confined to the mountain areas, resulting in a “ring” of higher species richness around the basin for some major lineages, including vascular plants, mammals, passeriform birds, and amphibians [28,29,30,31]. Some of these species are distributed continuously and may have a ring-shaped divergent pattern [32,33,34,35,36,37]. Many others have disjunct distributions ranging from much of the surrounding mountains to one area in the mountains, with most species distributions somewhere between these two extremes [38,39,40], which can potentially be explained by a dynamic history involving recurrent range fragmentation and contractions. In the western part of the region, the disjunct distribution of a variety of plant and animal species overlaps with the eastern edge of the Mountains of Southwest China hotspot [41,42,43,44,45,46,47,48], contributing to the biodiversity of the global hotspot and implying a shared long-term refugium [38, 39, 49,50,51,52,53]. Relevant studies that provided both species history inferences and multilocus estimates of genetic diversity distributions for such taxa are rare. They have focused on flowering plants, revealing local persistence through Quaternary climatic cycles and the possibility of extensive pollen-mediated gene flow while detecting, in hotspot edge areas, not necessarily unique or high genetic variation [39, 45, 54, 55]. Further case studies from other major evolutionary lineages would be beneficial.

Fig. 1
figure 1

Sampling sites of Leptobrachium boringii, the eastern part (red) of the Mountains of Southwest China biodiversity hotspot, and geographic pattern of expected heterozygosity (He) estimated from the 10,869-SNP data. Most known local distribution ranges are smaller than the sampling site symbols on the map. The others may cover wider areas enclosed with dotted lines. All the reported ranges were sampled. White indicates SNP data not collected. EE: the eastern edge of the hotspot. The map was created using DIVA-GIS version 7.5.0.0 (http://www.diva-gis.org/)

Full size image

The megophryid Emei mustache toad, Leptobrachium boringii (Liu, 1945), has a typical disjunct distribution in the mountains at the edges of the basin and hotspot (Fig. 1). This species occurs along the southern-half edge of the basin and in a small insular southern range, forming the northwestern range limit of the genus [56]. The toad inhabits mountain valleys covered by broad-leaved forest at elevations of 600 to 1700 m, breeds in permanent streams from February to March, and remains terrestrial outside the breeding season; additionally, the tadpoles usually take 3 years to reach metamorphosis [57, 58]. The “mustaches” are 10 to 16 maxillary keratinized nuptial spines used in male-male combat, and such structures occur only in L. boringii and close relatives among extant amphibians [58, 59]. Overcollecting for pet trade and habitat reduction are the main threats to this endangered (EN) species [60]. Previous studies with limited L. boringii sampling suggested pre-Quaternary divergences among major lineages [61] and a close relationship between the middle and eastern populations [62]. A further elaboration of the evolutionary processes and distribution of genetic diversity will contribute to the conservation of this species and may provide insights into the factors behind the similar disjunct distributions frequently observed in other frogs (e.g. Xenophrys omeimontis, Rhacophorus omeimontis, Babina daunchina, Bombina lichuanensis, [63, 64]).

In this study, we examine the distributional history and genetic diversity of the Emei mustache toad, focusing on two questions: (1) is the distribution highly dynamic over time, and what are the effects of range dynamics, and (2) does the eastern edge of the Mountains of Southwest China hotspot provide a long-term refugium harboring high genetic variation? We obtained sequence data for two mitochondrial loci and genome-wide single nucleotide polymorphism (SNP) data for samples from all of the corresponding known ranges. By combining time-calibrated gene tree, genetic diversity, genetic clustering, population splitting and mixing inference, and paleoclimatic suitability modeling, we inferred considerable range dynamics underlying the genetic and distribution patterns and a long-term refugium in the hotspot edge area, with the lowest genetic variability despite immigration. In addition, we demonstrated significant discordance between variability estimates from genome-wide SNPs and mitochondrial DNA sequences.

Results

Molecular data

The mitochondrial DNA dataset contained 1831 sites and 83 haplotypes after removing alignment gaps in the control region. Among the ingroup members, 282 sites were variable, and 227 sites were parsimony informative. Haplotype designation and GenBank accession numbers are given in Table S1. The genome-wide SNP main dataset contained 10,869 SNPs, from 7045 to 10,423 with a mean of 9305 SNPs (85.6%) for the ingroup and 6986 and 5680 SNPs for the two outgroup individuals. The 4-individual dataset with no missing loci included 16,888 SNPs. The Sequence Read Archive (SRA) accession numbers under BioProject PRJNA551927 are presented in Table S1. These three datasets are available on Dryad (https://0-doi-org.brum.beds.ac.uk/10.5061/dryad.t4b8gtj0c).

Genetic structure and divergence dates

For the BEAST analysis, a three-partition scheme was selected for the mitochondrial dataset: the protein-coding genes with the TrN + G model, tRNA genes with the K80 + I model, and control region with an HKY + I + G model. An adequate effective sample size of greater than 700 was achieved for each parameter [65], and the resulting topology was relatively well supported (Fig. 2). Haplotypes were rarely shared among sampling sites. Five major mitochondrial lineages that diverged at approximately 3.7 to 2.3 mya were detected, including four corresponding exclusively to regions South (sites 1 and 2), East (3–14), or West (21–28). In the other major lineage, haplotypes from the Middle region (sites 15–20) were nested within East clades, and sequential branching events from approximately 0.8 to 0.1 mya were inferred from East through Middle to West.

Fig. 2
figure 2

Time-calibrated mitochondrial gene tree of Leptobrachium boringii reconstructed with BEAST. Outgroups are not shown. Gray bars represent the 95% highest posterior densities for the age estimates. Nodes with Bayesian posterior probabilities ≥0.99, 0.95–0.98, and 0.80–0.94 are marked with black, gray, and open circles, respectively. Numbers appearing in haplotype names indicate sampling sites

Full size image

The first two principal coordinates obtained by PCoA jointly accounted for 71.9% of the observed genome-wide SNP variation among individuals (Fig. 3). Three discrete clusters corresponding to the South, East+Middle, and West regions were grouped. The latter two clusters were most closely related. Between them, samples from the most closely located sites 20 and 21 were also most closely related. The neighbor-joining tree (Fig. 3) could be summarized as (South,((East,Middle),West)), with individuals of each sampling site forming a unique cluster and an average bootstrap support of 97.2 for these clusters and their relationships. Sequential westward divergences were inferred among the Middle samples.

Fig. 3
figure 3

Principal coordinate plot and neighbor-joining tree (outgroups not shown) based on the SNP identity-by-state distance matrix for Leptobrachium boringii individuals. Numbers at nodes indicate bootstrap support values ≥90. Nodes with bootstrap support values of 100 and 97–99 are marked with closed and open circles, respectively. The names of the sampling sites correspond with those in Fig. 1. Outgroups are not shown

Full size image

Population relationships showing clear westward divergences and secondary contact were estimated by the TreeMix analysis of SNP data. Large residuals suggesting admixture between the Middle site 20 and the West sites, especially 21, were obtained when migrations were not allowed (Fig. 4). When migration edges were added, an edge with a weight of 0.476 from site 20 to the root of the West clade was found to most increase the likelihood. In addition, another three statistically significant edges (P < 0.05) involving some Middle, East, and South populations were added, with a resulting topology highly consistent with that of the neighbor-joining tree (Fig. 4).

Fig. 4
figure 4

Residual fit (top) from the maximum-likelihood population tree without migration and the admixture graph (bottom) with four migration edges (arrows) inferred by TreeMix. The names of the sampling sites correspond with those in Fig. 1. Outgroups are not shown. Large residuals indicated potential candidates for admixture events. The weight of the migration edge represents the fraction of ancestry derived from the edge

Full size image

Genetic diversity

Various levels of regional genetic variation and a decline in variability along the path of sequential westward divergence were inferred from the main SNP dataset (Fig. 1, Table 1). The expected heterozygosity and mean proportion of heterozygous loci estimates were similar and strongly correlated (Pearson’s r = 0.972, P < 0.001, n = 25), with the highest values for some South and East populations and lowest ones across the West region. Average identity-by-state distances of 0.032, 0.076, 0.075, and 0.058 were obtained for the West, Middle, East, and South regions, respectively. When compared to the outgroup species, again, much lower variation was estimated for the West representative. From the 4-individual SNP dataset without missing loci, heterozygous loci proportions of 0.0616, 0.1462, 0.1614, and 0.1152 were obtained for West representative bf71, East representative fjnB02, and the two outgroups ms07 and yb01, respectively.

Table 1 Genetic diversity estimates for each sampling site of Leptobrachium boringii. n: number of individuals; He: expected heterozygosity; \( \overline{H} \): mean proportion of heterozygous loci across individuals; h: number of haplotypes; Hd: haplotype diversity; π: nucleotide diversity
Full size table

A similar variability pattern was not obviously represented in the mitochondrial diversity measures (Table 1). The nucleotide diversity and haplotype diversity were not correlated with the SNP expected heterozygosity (P > 0.05), and only a weak correlation was detected between the heterozygosity and haplotype number (Pearson’s r = 0. 474, P = 0.017, n = 25). For example, for the five single-haplotype sampling sites, expected heterozygosities of 0.0300–0.0565 were obtained, covering a 56.1% range of the estimates. In 24 sampling sites, both mitochondrial and SNP data were obtained for at least two individuals. Similarly, when only these 74 individuals were considered, all the three mitochondrial diversity measures (results not shown) were not correlated with the SNP expected heterozygosity (P > 0.05).

Climatic suitability

The inferred regional climatic suitability of L. boringii in different periods was not equally stable (Fig. 5). The Maxent model provided reasonable discrimination, with an average AUC of 0.955 (SD = 0.015). The three variables that provided the greatest contributions to model development were the mean temperature of warmest quarter (37%), isothermality (25%), and precipitation of warmest quarter (13.9%). The remaining variables, including precipitation of driest month, annual precipitation, mean temperature of driest quarter, and temperature annual range, provided 8.9, 6.8, 4.4, and 4% contributions, respectively. Around West sites 21–26, the potential LIG, LGM, and MH distributions overlap extensively with the eastern edge of the Mountains of Southwest China hotspot. In contrast, at the southern edge of the Sichuan Basin, mountains harboring mainly the Middle sites were inferred to be generally unsuitable during the LIG and LGM periods. Compared with the distributions in the current interglacial (present and MH) period, a more restricted LIG projection was obtained.

Fig. 5
figure 5

Climatic suitability for Leptobrachium boringii determined by Maxent modeling. Circles are both occurrence records and sampling localities. The maps were created in the R environment using the package raster version 2.8 (https://cran.r-project.org/web/packages/raster/index.html)

Full size image

Discussion

Historical range dynamics and current disjunct distributions

In the mountains along the southern-half edge of the Sichuan Basin, our results indicate a Late Pliocene or early Pleistocene divergence between East and West L. boringii, much later westward dispersal establishing intermediate populations and causing secondary contact, and regionally varying climate stability. The East/West divergence was dated at 2.28 mya with a 95% highest posterior density of 3.45–1.28 mya (Fig. 2), overlapping with a previous estimate of 4.74–3.16 mya [61] and suggesting the long-term occupation of the two regions. Accordingly, a clear differentiation between the East and West samples was observed (Figs. 3 and 4). Westward dispersal occurring at approximately 0.8–0.1 mya was supported by the close relationship between the Middle and East samples, which agrees with the findings of Zheng et al. [62], and by the sequential westward divergences (Figs. 2, 3 and 4) accompanied by a significant loss of genome-wide SNP variation (Fig. 1, Table 1). The secondary contact is identified by the West mitochondrial DNA haplotypes nested within those of the western-most Middle site 20 (Fig. 2) and by a strong migration edge from site 20 to the West clade in the maximum-likelihood population tree inferred from the SNP data (Fig. 4). In addition to the pattern of secondary contact via dispersal after long-term isolation, the migration edges inferred within and between the East, Middle, and South regions suggest historical range changes. In line with this, local populations 4–8 are most closely related, and all exhibit a mitochondrial haplotype 4/5/6/7/8, implying recent admixture. It is likely that profound variations in regional climate stability, as exemplified by the modeled climatic suitability in the LIG, LGM, MH, and present-day periods (Fig. 5), have contributed to range dynamics and the current disjunct distribution of L. boringii. Nearly each part of the mountains surrounding the basin was inferred to be generally unsuitable in one or more periods. Meanwhile, distinct endemic genetic lineages such as those of sites 9, 13–16, 18, 20, and 26 (Fig. 2) support persistence in local refugia, which may be small given the complex mountain topography. We hence recommend caution in interpreting our climatic suitability results, and mainly use a general pattern implied by them together, recurrent range fragmentation and contraction, for hypothesis generating. Under suitable conditions, large range expansions might have not occurred in some refugial populations [8, 16]. This possibility and the recurrent range fragmentation and contraction caused by Quaternary climatic changes could lead to a largely disjunct distribution. Such a model provides a hypothesis for explaining scattered ranges and can be considered in future studies of the evolutionary histories of frogs and many other organisms isolated in the mountains surrounding the Sichuan Basin [38, 39, 48, 63, 64, 66, 67].

A substantially different modeled LIG distribution compared with those during the present interglacial period is not unique to the Emei mustache toad. Around the Sichuan Basin, such results have been commonly obtained for various animal [49, 50, 68,69,70,71] and plant groups [45, 55, 72, 73]. Because the LIG period at ~ 140–120 ka is relatively recent on an evolutionary scale, it is unlikely that the corresponding effects have been erased. In fact, a few studies have highlighted the significant roles of the LIG climate in shaping species histories [49, 69, 74, 75]. In theory, under global climatic cycles, profound climatic differences between the last and present interglacial periods should be limited to certain regions, which vary with the climatic variables of interest. An awareness of this relation may be beneficial, especially for examining the evolutionary history of narrow-ranged species and understanding regional biodiversity patterns.

The basal phylogenetic position of the southern population (Figs. 2 and 4) provides a clue to the early distributional history of the species. In the present distribution of L. boringii, the South region (Fig. 1) is most close to the distribution or proposed ancestral range of the close relatives of the species [15, 56, 62, 76,77,78]. It may be not far away from the ancestral range of the species and, as suggested by the high and unique genetic variation estimates (Fig. 3, Table 1), has probably served as a long-term refugium. One plausible scenario is that the eastern range was established through northward expansion no later than early Pleistocene, given that the West clade with a Late Pliocene to early Pleistocene origin was nested within eastern lineages according to our time-calibrated tree (Fig. 2). Then, the eastern population might have expanded to the western range by a route along mountains at the southern edge of the basin, which were reoccupied through the aforementioned, much latter westward dispersal.

Spatial heterogeneity of the factors that influence high biodiversity within a hotspot

This study detects both long-term refugial and young immigrant (site 20) populations in the eastern edge of the Mountains of Southwest China hotspot. The molecular results suggest a Late Pliocene to early Pleistocene origin for the West clade endemic to the hotspot edge, in which the area around sites 21–24 and 26 has consistently been modeled to be suitable for L. boringii in the LIG, LGM, MH, and current periods. Western marginal refugia and areas of long-term persistence within or overlapping with this hotspot edge have been proposed for some mountain amphibian [37, 79], mammal [47], insect [50], and flowering plant species [38, 39, 52, 53, 80,81,82]. Additionally, gene flows from the east via mountains along the edges of the Sichuan Basin have recently been reported for different species [37, 48, 83]. Together with these findings, our results imply that the preservation of eastern taxa, facilitated by the continuous mountain topography surrounding the basin, contributes significantly to regional biodiversity.

In this hotspot, except for the southwestern finger, which overlaps the Nujiang Valley, the eastern edge is the only area that provides low-elevation habitat, say, below 1000 m. Most areas of the hotspot exceed 3000 m due to the rapid uplift of the Hengduan Mountain region mainly between the late Miocene and Late Pliocene [24, 84,85,86,87], which has been identified as a major driver of the diversification of hotspot flora by creating conditions favoring in situ speciation [88]. These areas are generally much higher than the mainland to the east and are not suitable for low-elevation species. Indeed, in addition to the aforementioned taxa, the eastern edge of the hotspot forms a western range limit for many other species [45, 54, 64, 89,90,91]; this is not surprising given that the biodiversity hotspots in temperate realms usually host a small proportion of endemics and numerous immigrant species, especially for mammals and birds [92, 93]. However, the pattern suggests that a certain contributor to biodiversity, preserving and receiving low-elevation elements in this case, can be significant for only a particular part of a hotspot. In noting such heterogeneity, we do not overlook the fact that the eastern edge exhibits highly complex geological structures that facilitate diversification [71] and harbors endemics at different elevations [94, 95]. Hotspots are critical for global biodiversity conservation [92]. Further studies on the heterogeneous nature of each hotspot may aid in conservation planning.

Genome-wide SNP versus mitochondrial DNA diversity estimates

In addition to implying that long-term occupation does not necessarily correlate with large genetic variation at the western edge of the Sichuan Basin (Figs. 1, 2 and 5), our results provide an instance of discordant variation estimates using genome-wide SNPs and mitochondrial DNA sequences. Among the mitochondrial results, only the haplotype numbers are weakly correlated with the SNP expected heterozygosities (Table 1). Relatively reliable estimates of heterozygosity can be obtained using a large number of SNP loci from even a small number of individuals [20]. Compared with such estimates, mitochondrial DNA diversity is more susceptible to genetic drift and selective sweeps due to the smaller effective population size and lack of recombination. Moreover, for the same reason and because of possible sex-biased dispersal, population admixture may not have left proportional immigrant mitochondrial alleles. This uncertainty can also lead to discordant estimates from the two types of data. A previous study of this species has reported genetic parameters consistent with female-biased dispersal [96]. In the West region, site 25 was unique for having mitochondrial haplotype samples of both local and immigrant origins (Fig. 2) and an estimated nucleotide diversity of 0.0111. This value was at least five times larger than those of the other West sites and was higher than estimates for most sites in the other regions. However, the SNP heterozygosity estimates for West sites were similar and much lower than those in all the other regions. Given the high genetic homogeneity revealed by a substantially smaller average identity-by-state distance (0.032) than that in other regions (0.058–0.076), we suggest considering a bottleneck in addition to evolutionary stasis and other scenarios in an effort to explain the low variability found in western L. boringii.

Different RADseq and mitochondrial variations associated with admixture have been reported in other taxa [97,98,99]. The discordance detected in the present work is relevant to the advantage of using a large number of loci to estimate levels of genetic variation [20, 21]; notably, different estimates from genome-wide loci and mitochondrial DNA can be obtained for populations that experienced admixture. An awareness of this possibility may facilitate the integration of results from different sources of data.

Conclusions

A scenario of dispersal after long-term isolation can explain much of the spatial distribution of genetic variation in the Emei mustache toad, yet long-term persistence and secondary contact do not guarantee high variability. For many other taxa isolated in the mountains surrounding the Sichuan Basin, a model that combines recurrent contractions caused by Quaternary climatic changes and some failed expansions under suitable conditions for explaining scattered ranges can be considered when examining their evolutionary histories. More generally, this study provides insights into the heterogeneous nature of biodiversity hotspots and the discordant variability estimates using genome-wide and mitochondrial data.

Methods

Sampling

A total of 301 individuals, including 289 tadpoles, from 28 sampling sites representing all the known ranges of Leptobrachium boringii were used as the ingroup. The congeners L. liui and L. leishanense, with two tadpole samples each, were used as outgroups according to the current phylogenetic hypotheses [62, 76,77,78]. Sequences of two mitochondrial DNA fragments were obtained for all samples, and genome-wide SNPs were obtained for a representative subset of 77 individuals. Detailed sampling information is provided in Table S1.

Mitochondrial DNA data

The 745-bp nad1 fragment consisted of partial tRNA-Leu (16 bp) and NADH dehydrogenase subunit 1 genes; the D-loop fragment of ~ 1080–1350 bp in length contained regions of the cytochrome b gene (78 bp) and a control region with the tRNA-Thr, tRNA-Pro, and tRNA-Trp genes between them. For the nad1 fragment, the primers Leu-1 L and Leu-5 L from Zheng et al. [62] and ND1-6H from Zheng et al. [15] were used in PCR and Sanger sequencing (ABI 3730). Sequences of this fragment for 19 individuals were also obtained from these two previous studies. The primers for the D-loop fragment were newly designed in this study: Vcb-1 L (5′-ATCggCggAC AACCTgTCgA AgACC-3′), Vcb-3 L (5′-CAACCTgTCg AAgACCCCTA CgT-3′), Vcr-2H (5′-gTggCTgATC CACCggAAgg TAAg-3′), and Vcr-4H (5′-gCTgATCCAC CggAAggTAA gATC-3′). Sequences of the light-strand encoded tRNA-Pro gene were converted into complementary strand sequences. Alignment was conducted with ClustalX version 2.1 [100] and visually checked.

Mitochondrial diversity and time-calibrated gene tree

The nad1 and D-loop fragments were concatenated for analysis. The nucleotide diversity and haplotype diversity were estimated for each local population using DnaSP version 6.10.01 [101]. To examine the population relationships and the dates of divergence events, a time-calibrated tree was reconstructed in a Bayesian framework with BEAST version 2.5.1 [102]. The best-fitting combination of nucleotide partitions and substitution models was determined with PartitionFinder version 2.1.1 [103] using the Bayesian information criterion. To avoid potential over-parameterization of the intraspecific analysis, only three data blocks were defined for evaluation: the protein-coding genes, tRNA genes, and control region. The protein-coding gene block contained mainly the nad1 gene sequence. Because suitable fossil calibrations were not available, an evolutionary rate of 0.69 ± 0.3% divergence per million years per lineage commonly used for the anuran nad1 gene was used to calibrate the tree [15, 104,105,106]. This approach required the protein-coding gene block to be constrained as one partition in analysis. The uncorrelated lognormal model was used in the relaxed clock method. The calibrated Yule model was specified for the tree prior. The uncorrelated lognormal relaxed clock mean (ucldMean) was constrained between 0.0039 and 0.0099 substitutions per site per million years, and a normal distribution was applied (mean 0.0069, sigma 0.0015). Four independent runs of 20 million generations each were performed to avoid local optima, and a sampling frequency of 2000 generations was applied. The performance of the runs was evaluated using Tracer version 1.6 [107] to ensure that the chains were converging and mixing sufficiently, and the final 90% of samples for each of the four runs were combined to produce a maximum clade credibility tree.

Restriction site-associated DNA sequencing and SNP data

Genome-wide SNP data from RADseq loci were obtained for 75 ingroup representatives from 25 sampling sites (one to seven individuals per site, mostly three) and two outgroup individuals using the SLAF-seq [108] approach. Genomic DNA (~ 0.15 μg per sample) was incubated at 37 °C with the MseI restriction enzyme, T4 DNA ligase, ATP, and MseI adaptor containing barcode information. The restriction ligation was heat inactivated at 65 °C. The product was digested with the restriction enzyme HaeII at 37 °C and subsequently purified using Agencourt AMPure XP beads (Beckman Coulter). PCR was performed using the purified product and Phusion Master Mix (New England Biolabs) with universal and index primers. The PCR products were purified and pooled and then electrophoresed on 2% agarose gel. Fragments of 425–450 bp in size (with indexes and adaptors) were isolated, purified and diluted for sequencing using the PE150 Illumina Hiseq 2000 platform.

Raw reads were checked for quality using FastQC version 0.11.7 [109]. Data cleaning was conducted using the process_radtags program in Stacks version 2.0b [110], discarding reads with an uncalled base or an average Phred33 score of < 10 within a 22-bp sliding window. A de novo assembly of loci was conducted using the denovo_map.pl pipeline in Stacks. In a preliminary analysis of a representative subset of 32 individuals, the minimum number of reads required to form an allele, number of mismatches allowed between alleles, and number of mismatches allowed between loci during the construction of the catalog in Stacks were determined as 3, 4, and 5, respectively. These values maximized the number of polymorphic loci found in 80% of the individuals [111]. Through this pipeline, biallelic SNPs present in at least 80% of individuals with a minor allele frequency greater than 1%, one random SNP per locus, were obtained for further filtering using VCFtools version 0.1.17 [112]. The final main dataset was composed of SNPs with a mean read depth across individuals greater than 6 and below 50 to avoid the potential effects of sequencing errors and repetitive regions. In addition, keeping all settings equal except for allowing no missing data, a 4-individual dataset was generated to compare the genetic diversity of two ingroup representatives and two outgroup individuals.

Genetic diversity and structure analyses of SNPs

The expected heterozygosity for each sampling site was estimated using Arlequin version 3.5.2.2 [113]. In addition, the heterozygosity for each individual was measured by the proportion of heterozygous loci using the --het option in VCFtools and then averaged across individuals for each sampling site. Using SPSS version 12.0, various diversity estimates were first checked for normality with the Kolmogorov-Smirnov test, and then a Pearson or Spearman correlation test was performed for comparison.

An identity-by-state distance matrix between individuals was constructed in PLINK version 1.9 [114] and used for principal coordinates analysis (PCoA) and neighbor-joining tree building. PCoA was conducted using GenAlEx version 6.502 [115]. A neighbor-joining tree was constructed using the Neighbor program in PHYLIP version 3.695 [116], with 500 bootstrap replicates to evaluate the nodal support. A Perl script VCFbootstrap.pl (https://0-doi-org.brum.beds.ac.uk/10.5061/dryad.t4b8gtj0c) was written to generate the 500 bootstrap SNP datasets in VCF format. At the population level, a maximum-likelihood inference of splitting and mixing from the allele frequency was conducted using TreeMix version 1.13 [117] to assess the diversification and potential secondary contacts. The allele frequency at the sampling site level was generated by PLINK. Models allowing zero to five migration events were explored in separate TreeMix analyses, and the results were compared.

Climatic suitability

As an approach to assess range stability, the last interglacial (LIG), last glacial maximum (LGM), and mid-Holocene (MH) potential distributions were modeled for comparison. Climate reconstructions for these periods are frequently used examples of Quaternary climate fluctuations. The analyses were conducted at a spatial extent of 21.5–34°N and 100–114°E, extending not substantially from the Sichuan Basin and surrounding mountains and the range of L. boringii.

Ecological niche modeling was conducted using Maxent version 3.4.1 [118, 119]. Nineteen current bioclimatic variables with a resolution of 2.5 arc-minutes were obtained from the WorldClim database [120]. The occurrence data containing 28 records, with no more than one occurrence per grid cell, were derived from our field surveys (https://0-doi-org.brum.beds.ac.uk/10.5061/dryad.t4b8gtj0c). For highly correlated temperature/precipitation variable pairs (Pearson, |r| > 0.8), the variables with smaller contributions according to jackknife analysis of variable importance to model development were excluded [118, 121, 122]. Seven variables were retained: isothermality, temperature annual range, mean temperature of driest quarter, mean temperature of warmest quarter, annual precipitation, precipitation of driest month, and precipitation of warmest quarter. In a Maxent analysis, we ran models with 10 bootstrap replicates by randomly assigning occurrence records to the training (75%) and testing (25%) data sets. The replicates were used as proxy models to develop consensus-based ensemble forecasts [123], and the mean suitability was output in logistic format [124]. Model performance was assessed using the area under the receiver operating curve (AUC).

The niche model was projected onto paleoclimate reconstructions in WorldClim with the same resolution of 2.5 arc-minutes. The resolution of the LIG (~ 140–120 thousand years ago; ka) layers was changed from 30 arc-seconds to 2.5 arc-minutes through averaging [125]. For LGM (22 ka) and MH (6 ka), taking the uncertainty stemming from different global climate models into consideration, the commonly used reconstructions by the Community Climate System Model (CCSM4, [126]) and the Model for Interdisciplinary Research on Climate (MIROC-ESM, [127]) were adopted.

Availability of data and materials

Mitochondrial DNA sequences: GenBank accessions KP054478–KP054954 and MK164282–MK164395. RADseq raw reads: Sequence Read Archive accessions SRR9645625–SRR9645701 under BioProject PRJNA551927. The Perl script VCFbootstrap.pl, occurrence data for ecological niche modeling, mitochondrial DNA dataset, main SNP dataset, and 4-individual SNP dataset are available on Dryad (https://0-doi-org.brum.beds.ac.uk/10.5061/dryad.t4b8gtj0c).

Abbreviations

ATP:

Adenosine triphosphate

AUC:

Area under the receiver operating curve

CCSM:

Community Climate System Model

LGM:

Last glacial maximum

LIG:

Last interglacial

MH:

Mid-Holocene

MIROC:

Model for Interdisciplinary Research on Climate

mtDNA:

Mitochondrial deoxyribonucleic acid

PCoA:

Principal coordinates analysis

PCR:

Polymerase chain reaction

RADseq:

Restriction-site associated deoxyribonucleic acid

SLAF-seq:

Specific-locus amplified fragment sequencing

SNP:

Single nucleotide polymorphism

VCF:

Variant call format

References

  1. Moritz C. Strategies to protect biological diversity and the evolutionary processes that sustain it. Syst Biol. 2002;51:238–54.

    PubMed  Google Scholar 

  2. Hewitt GM. The structure of biodiversity–insights from molecular phylogeography. Front Zool. 2004;1:4.

  3. Palsbøll PJ, Berube M, Allendorf FW. Identification of management units using population genetic data. Trends Ecol Evol. 2007;22:11–6.

    PubMed  Google Scholar 

  4. Hewitt GM. Genetic consequences of climatic oscillations in the Quaternary. Phil Trans R Soc B Biol Sci. 2004;359:183–95.

  5. Hofreiter M, Stewart J. Ecological change, range fluctuations and population dynamics during the Pleistocene. Curr Biol. 2009;19:R584–94.

    CAS  PubMed  Google Scholar 

  6. Stewart JR, Lister AM, Barnes I, Dalén L. Refugia revisited: individualistic responses of species in space and time. Proc R Soc B Biol Sci. 2010;277:661–71.

  7. Hewitt G. The genetic legacy of the Quaternary ice ages. Nature. 2000;405:907–13.

  8. Bennett KD, Provan J. What do we mean by ‘refugia’? Quat Sci Rev. 2008;27:2449–55.

    Google Scholar 

  9. Petit RJ, Aguinagalde I, de Beaulieu JL, Bittkau C, Brewer S, Cheddadi R, et al. Glacial refugia: hotspots but not melting pots of genetic diversity. Science. 2003;300:1563–5.

    CAS  PubMed  Google Scholar 

  10. Teixeira J, Gonçalves H, Ferrand N, García-París M, Recuero E. Mitochondrial phylogeography of the Iberian endemic frog Rana iberica, with implications for its conservation. Current Zoology. 2018;64:755–64.

    PubMed  PubMed Central  Google Scholar 

  11. Carnaval AC, Hickerson MJ, Haddad CF, Rodrigues MT, Moritz C. Stability predicts genetic diversity in the Brazilian Atlantic forest hotspot. Science. 2009;323:785–9.

    CAS  PubMed  Google Scholar 

  12. Yannic G, Pellissier L, Ortego J, Lecomte N, Couturier S, Cuyler C, et al. Genetic diversity in caribou linked to past and future climate change. Nat Clim Chang. 2014;4:132–7.

    Google Scholar 

  13. Hu J, Huang Y, Jiang J, Guisan A. Genetic diversity in frogs linked to past and future climate change on the roof of the world. J Anim Ecol. 2019;88:953–63.

    PubMed  Google Scholar 

  14. Canestrelli D, Cimmaruta R, Costantini V, Nascetti G. Genetic diversity and phylogeography of the Apennine yellow-bellied toad Bombina pachypus, with implications for conservation. Mol Ecol. 2006;15:3741–54.

    CAS  PubMed  Google Scholar 

  15. Zheng Y, Hu J, Zeng X. Examining the interglacial high-elevation refugia scenario in East Asian subtropical mountain systems with the frog species Leptobrachium liui. Ecol Evol. 2018;8:9326–40.

  16. Gómez A, Lunt DH. Refugia within refugia: patterns of phylogeographic concordance in the Iberian Peninsula. In: Weiss S, Ferrand N, editors. Phylogeography of southern European refugia. Dordrecht: Springer; 2007. p. 155–88.

    Google Scholar 

  17. Abellán P, Svenning JC. Refugia within refugia–patterns in endemism and genetic divergence are linked to Late Quaternary climate stability in the Iberian Peninsula. Biol J Linn Soc. 2014;113:13–28.

    Google Scholar 

  18. Gavin DG, Fitzpatrick MC, Gugger PF, Heath KD, Rodríguez-Sánchez F, Dobrowski SZ, et al. Climate refugia: joint inference from fossil records, species distribution models and phylogeography. New Phytol. 2014;204:37–54.

    PubMed  Google Scholar 

  19. Roberts DR, Hamann A. Glacial refugia and modern genetic diversity of 22 western North American tree species. Proc R Soc B Biol Sci. 2015;282:20142903.

  20. Nei M. Estimation of average heterozygosity and genetic distance from a small number of individuals. Genetics. 1978;89:583–90.

    CAS  PubMed  PubMed Central  Google Scholar 

  21. Garrick RC, Bonatelli IA, Hyseni C, Morales A, Pelletier TA, Perez MF, Rice E, Satler JD, Symula RE, Thome MTC, Carstens BC. The evolution of phylogeographic data sets. Mol Ecol. 2015;24:1164–71.

    CAS  PubMed  Google Scholar 

  22. Clark MK, Schoenbohm LM, Royden LH, Whipple KX, Burchfiel BC, Zhang X, Tang W, Wang E, Chen L. Surface uplift, tectonics, and erosion of eastern Tibet from large-scale drainage patterns. Tectonics. 2004;23:TC1006.

    Google Scholar 

  23. Shen C, Mei L, Xu Z, Tang J. Architecture and tectonic evolution of composite basin-mountain system in Sichuan basin and its adjacent areas. Geotecton Metallog. 2007;31:288–99.

    Google Scholar 

  24. Sun BN, Wu JY, Liu YS, Ding ST, Li XC, Xie SP, Yan DF, Lin ZC. Reconstructing Neogene vegetation and climates to infer tectonic uplift in western Yunnan, China. Palaeogeogr Palaeoclimatol Palaeoecol. 2011;304:328–36.

    Google Scholar 

  25. Liu S, Deng B, Li Z, Sun W. Architecture of basin-mountain systems and their influences on gas distribution: a case study from the Sichuan basin, South China. J Asian Earth Sci. 2012;47:204–15.

    Google Scholar 

  26. Li JJ, Shu Q, Zhou SZ, Zhao ZJ, Zhang JM. Review and prospects of Quaternary glaciation research in China. J Glaciol Geocryol. 2004;26:235–43.

  27. Owen LA, Dortch JM. Nature and timing of Quaternary glaciation in the Himalayan–Tibetan orogen. Quat Sci Rev. 2014;88:14–54.

  28. Mutke J, Barthlott W. Patterns of vascular plant diversity at continental to global scales. Biologiske Skrifter. 2005;55:521–31.

    Google Scholar 

  29. Roselaar CS, Sluys R, Aliabadian M, Mekenkamp PGM. Geographic patterns in the distribution of Palearctic songbirds. J Ornithol. 2007;148:271–80.

    Google Scholar 

  30. López-Pujol J, Zhang FM, Sun HQ, Ying TS, Ge S. Centres of plant endemism in China: places for survival or for speciation? J Biogeogr. 2011;38:1267–80.

    Google Scholar 

  31. Jenkins CN, Pimm SL, Joppa LN. Global patterns of terrestrial vertebrate diversity and conservation. Proc Natl Acad Sci. 2013;110:E2602–10.

    CAS  PubMed  PubMed Central  Google Scholar 

  32. Monahan WB, Pereira RJ, Wake DB. Ring distributions leading to species formation: a global topographic analysis of geographic barriers associated with ring species. BMC Biol. 2012;10:20.

    PubMed  PubMed Central  Google Scholar 

  33. He K, Jiang X. Sky islands of Southwest China. I: an overview of phylogeographic patterns. Chin Sci Bull. 2014;59:585–97.

    Google Scholar 

  34. Pan T, Wang H, Hu C, Sun Z, Zhu X, Meng T, Meng X, Zhang B. Species delimitation in the genus Moschus (Ruminantia: Moschidae) and its high-plateau origin. PLoS One. 2015;10:e0134183.

    PubMed  PubMed Central  Google Scholar 

  35. Kuchta SR, Wake DB. Wherefore and whither the ring species? Copeia. 2016;104:189–201.

    Google Scholar 

  36. Yang LQ, Hu HY, Xie C, Lai SP, Yang M, He XJ, Zhou SD. Molecular phylogeny, biogeography and ecological niche modelling of Cardiocrinum (Liliaceae): insights into the evolutionary history of endemic genera distributed across the Sino-Japanese floristic region. Ann Bot. 2017;119:59–72.

    PubMed  Google Scholar 

  37. Qiao L, Wen G, Qi Y, Lu B, Hu J, Song Z, Fu J. Evolutionary melting pots and reproductive isolation: a ring-shaped diversification of an odorous frog (Odorrana margaratea) around the Sichuan Basin. Mol Ecol. 2018;27:4888–900.

    PubMed  Google Scholar 

  38. Xie XF, Yan HF, Wang FY, Ge XJ, Hu CM, Hao G. Chloroplast DNA phylogeography of Primula ovalifolia in central and adjacent southwestern China: past gradual expansion and geographical isolation. J Syst Evol. 2012;50:284–94.

    Google Scholar 

  39. Ma Q, Du YJ, Chen N, Zhang LY, Li JH, Fu CX. Phylogeography of Davidia involucrata (Davidiaceae) inferred from cpDNA haplotypes and nSSR data. Syst Bot. 2015;40:796–810.

    Google Scholar 

  40. Li S, Xu N, Lv J, Jiang J, Wei G, Wang B. A new species of the odorous frog genus Odorrana (Amphibia, Anura, Ranidae) from southwestern China. PeerJ. 2018;6:e5695.

    PubMed  PubMed Central  Google Scholar 

  41. Myers N, Mittermeier RA, Mittermeier CG, da Fonseca GAB, Kent J. Biodiversity hotspots for conservation priorities. Nature. 2000;403:853–8.

    CAS  PubMed  Google Scholar 

  42. Mittermeier RA, Gil PR, Hoffmann M, Pilgrim J, Brooks T, Mittermeier CG, Lamoreux J, da Fonseca GAB. Hotspots revisited: earth’s biologically richest and most endangered terrestrial ecoregions. Mexico City: CEMEX; 2004.

    Google Scholar 

  43. Qiu YX, Guan BC, Fu CX, Comes HP. Did glacials and/or interglacials promote allopatric incipient speciation in East Asian temperate plants? Phylogeographic and coalescent analyses on refugial isolation and divergence in Dysosma versipellis. Mol Phylogenet Evol. 2009;51:281–93.

  44. Zheng Y, Fu J, Li S. Toward understanding the distribution of Laurasian frogs: a test of Savage’s biogeographical hypothesis using the genus Bombina. Mol Phylogenet Evol. 2009;52:70–83.

    CAS  PubMed  Google Scholar 

  45. Bai WN, Wang WT, Zhang DY. Contrasts between the phylogeographic patterns of chloroplast and nuclear DNA highlight a role for pollen-mediated gene flow in preventing population divergence in an East Asian temperate tree. Mol Phylogenet Evol. 2014;81:37–48.

  46. Wei X, Sork VL, Meng H, Jiang M. Genetic evidence for central-marginal hypothesis in a Cenozoic relict tree species across its distribution in China. J Biogeogr. 2016;43:2173–85.

    Google Scholar 

  47. Lv X, Cheng J, Meng Y, Chang Y, Xia L, Wen Z, Ge D, Liu S, Yang Q. Disjunct distribution and distinct intraspecific diversification of Eothenomys melanogaster in South China. BMC Evol Biol. 2018;18:50.

    PubMed  PubMed Central  Google Scholar 

  48. Ye Z, Yuan J, Li M, Damgaard J, Chen P, Zheng C, Yu H, Fu S, Bu W. Geological effects influence population genetic connectivity more than Pleistocene glaciations in the water strider Metrocoris sichuanensis (Insecta: Hemiptera: Gerridae). J Biogeogr. 2018;45:690–701.

    Google Scholar 

  49. Wang W, McKay BD, Dai C, Zhao N, Zhang R, Qu Y, Song G, Li SH, Liang W, Yang X, Pasquet E, Lei F. Glacial expansion and diversification of an East Asian montane bird, the green-backed tit (Parus monticolus). J Biogeogr. 2013;40:1156–69.

  50. Ye Z, Zhu G, Chen P, Zhang D, Bu W. Molecular data and ecological niche modelling reveal the Pleistocene history of a semi-aquatic bug (Microvelia douglasi douglasi) in East Asia. Mol Ecol. 2014;23:3080–96.

    CAS  PubMed  Google Scholar 

  51. Guan BC, Chen W, Gong X, Wu T, Cai QY, Liu YZ, Ge G. Landscape connectivity of Cercidiphyllum japonicum, an endangered species and its implications for conservation. Ecological Informatics. 2016;33:51–6.

    Google Scholar 

  52. Bai G, Zhou T, Zhang X, Chen X, Yang J, Li Z, Zhao G. Genetic differentiation and population genetic structure of the Chinese endemic Dipteronia Oliv. revealed by cpDNA and AFLP data. Forests. 2017;8:424.

  53. Wang YH, Comes HP, Cao YN, Guo R, Mao YR, Qiu YX. Quaternary climate change drives allo-peripatric speciation and refugial divergence in the Dysosma versipellis-pleiantha complex from different forest types in China. Sci Rep. 2017;7:40261.

  54. Ying LX, Zhang TT, Chiu CA, Chen TY, Luo SJ, Chen XY, Shen ZH. The phylogeography of Fagus hayatae (Fagaceae): genetic isolation among populations. Ecol Evol. 2016;6:2805–16.

    PubMed  PubMed Central  Google Scholar 

  55. Zhang YH, Wang IJ, Comes HP, Peng H, Qiu YX. Contributions of historical and contemporary geographic and environmental factors to phylogeographic structure in a Tertiary relict species, Emmenopterys henryi (Rubiaceae). Sci Rep. 2016;6:24041.

  56. Frost DR. Amphibian species of the world: an online reference. Version 6.0. American Museum of Natural History. 2019. http://research.amnh.org/herpetology/amphibia/index.html. Accessed 28 May 2019.

  57. Wu G, Yang W. Studies on genus Vibrissaphora (Amphibia: Pelobatidae) 2. Some ecological notes of vibrissaphorids. Acta Herpetologica Sinica (Old Series). 1981;5:77–80.

    Google Scholar 

  58. Zheng Y, Rao D, Murphy RW, Zeng X. Reproductive behavior and underwater calls in the Emei mustache toad, Leptobrachium boringii. Asian Herpetological Research. 2011;2:199–215.

    Google Scholar 

  59. Hudson CM, He X, Fu J. Keratinized nuptial spines are used for male combat in the Emei moustache toad (Leptobrachium boringii). Asian Herpetological Research. 2011;2:142–8.

    Google Scholar 

  60. Fei L, Wu G. Leptobrachium boringii. The IUCN red list of threatened species 2004: e.T57625A11665713. 2004. https://0-doi-org.brum.beds.ac.uk/10.2305/IUCN.UK.2004.RLTS.T57625A11665713.en. Accessed 29 May 2019.

  61. Li J, Wei S, Hu M, Luo Z, Zhao M, Wu H. Reflection of paleoclimate oscillations and tectonic events in the phylogeography of moustache toads in southern China. J Zool. 2018;305:17–26.

    Google Scholar 

  62. Zheng Y, Li S, Fu J. A phylogenetic analysis of the frog genera Vibrissaphora and Leptobrachium, and the correlated evolution of nuptial spine and reversed sexual size dimorphism. Mol Phylogenet Evol. 2008;46:695–707.

    CAS  PubMed  Google Scholar 

  63. Fei L, Hu S, Ye C, Huang Y. Fauna Sinica. Amphibia. Volume 2. Anura. Beijing: Science Press; 2009.

  64. Fei L, Hu S, Ye C, Huang Y. Fauna Sinica. Amphibia. Volume 3. Anura Ranidae. Beijing: Science Press; 2009.

  65. Drummond AJ, Ho SYW, Phillips MJ, Rambaut A. Relaxed phylogenetics and dating with confidence. PLoS Biol. 2006;4:e88.

    PubMed  PubMed Central  Google Scholar 

  66. Gu L, Liu Y, Que P, Zhang Z. Quaternary climate and environmental changes have shaped genetic differentiation in a Chinese pheasant endemic to the eastern margin of the Qinghai-Tibetan Plateau. Mol Phylogenet Evol. 2013;67:129–39.

  67. Ge D, Lu L, Xia L, Chang Y, Wen Z, Lv X, Du Y, Liu Q, Yang Q. An endemic rat species complex is evidence of moderate environmental changes in the terrestrial biodiversity centre of China through the late Quaternary. Sci Rep. 2017;7:46127.

  68. He K, Hu NQ, Chen X, Li JT, Jiang XL. Interglacial refugia preserved high genetic diversity of the Chinese mole shrew in the mountains of Southwest China. Heredity. 2016;116:23–32.

    CAS  PubMed  Google Scholar 

  69. Dong F, Hung CM, Li XL, Gao JY, Zhang Q, Wu F, Lei FM, Li SH, Yang XJ. Ice age unfrozen: severe effect of the last interglacial, not glacial, climate change on East Asian avifauna. BMC Evol Biol. 2017;17:244.

  70. Wang P, Yao H, Gilbert KJ, Lu Q, Hao Y, Zhang Z, Wang N. Glaciation-based isolation contributed to speciation in a Palearctic alpine biodiversity hotspot: evidence from endemic species. Mol Phylogenet Evol. 2018;129:315–24.

    PubMed  Google Scholar 

  71. Zhao M, Chang Y, Kimball RT, Zhao J, Lei F, Qu Y. Pleistocene glaciation explains the disjunct distribution of the chestnut-vented nuthatch (Aves, Sittidae). Zool Scr. 2019;48:33–45.

    Google Scholar 

  72. Fan L, Zheng H, Milne RI, Zhang L, Mao K. Strong population bottleneck and repeated demographic expansions of Populus adenopoda (Salicaceae) in subtropical China. Ann Bot. 2018;121:665–79.

    PubMed  PubMed Central  Google Scholar 

  73. Niu YT, Ye JF, Zhang JL, Wan JZ, Yang T, Wei XX, Lu LM, Li JH, Chen ZD. Long-distance dispersal of postglacial contraction? Insights into disjunction between Himalaya-Hengduan Mountains and Taiwan in a cold-adapted herbaceous genus, Triplostegia. Ecol Evol. 2018;8:1131–46.

    PubMed  Google Scholar 

  74. Gür H. The effects of the Late Quaternary glacial–interglacial cycles on Anatolian ground squirrels: range expansion during the glacial periods? Biol J Linn Soc. 2013;109:19–32.

    Google Scholar 

  75. Perktaş U, Gür H, Sağlam İK, Quintero E. Climate-driven range shifts and demographic events over the history of Kruper's nuthatch Sitta krueperi. Bird Study. 2015;62:14–28.

    Google Scholar 

  76. Matsui M, Hamidy A, Murphy RW, Khonsue W, Yambun P, Shimada T, et al. Phylogenetic relationships of megophryid frogs of the genus Leptobrachium (Amphibia, Anura) as revealed by mtDNA gene sequences. Mol Phylogenet Evol. 2010;56:259–72.

    PubMed  Google Scholar 

  77. Pyron RA. Biogeographic analysis reveals ancient continental vicariance and recent oceanic dispersal in amphibians. Syst Biol. 2014;63:779–97.

    PubMed  Google Scholar 

  78. Jetz W, Pyron RA. The interplay of past diversification and evolutionary isolation with present imperilment across the amphibian tree of life. Nature Ecol Evol. 2018;2:850–8.

    Google Scholar 

  79. Li J, Zhao M, Wei S, Luo Z, Wu H. Geologic events coupled with Pleistocene climatic oscillations drove genetic variation of Omei treefrog (Rhacophorus omeimontis) in southern China. BMC Evol Biol. 2015;15:289.

    CAS  PubMed  PubMed Central  Google Scholar 

  80. Qi XS, Chen C, Comes HP, Sakaguchi S, Liu YH, Tanaka N, Sakio H, Qiu YX. Molecular data and ecological niche modelling reveal a highly dynamic evolutionary history of the East Asian Tertiary relict Cercidiphyllum (Cercidiphyllaceae). New Phytol. 2012;196:617–30.

  81. Wang JF, Gong X, Chiang YC, Kuroda C. Phylogenetic patterns and disjunct distribution in Ligularia hodgsonii Hook. (Asteraceae). J Biogeogr. 2013;40:1741–54.

  82. Cao YN, Comes HP, Sakaguchi S, Chen LY, Qiu YX. Evolution of East Asia’s Arcto-Tertiary relict Euptelea (Eupteleaceae) shaped by Late Neogene vicariance and Quaternary climate change. BMC Evol Biol. 2016;16:66.

  83. Xue J, Zhang H, Ning X, Bu W, Yu X. Evolutionary history of a beautiful damselfly, Matrona basilaris, revealed by phylogeographic analyses: the first study of an odonate species in mainland China. Heredity. 2019;122:570–81.

    PubMed  Google Scholar 

  84. Kirby E, Reiners PW, Krol MA, Whipple KX, Hodges KV, Farley KA, Tang W, Chen Z. Late Cenozoic evolution of the eastern margin of the Tibetan Plateau: inferences from 40Ar/39Ar and (U-Th)/He thermochronology. Tectonics. 2002;21(1):1.

  85. Clark MK, House MA, Royden LH, Whipple KX, Burchfiel BC, Zhang X, Tang W. Late Cenozoic uplift of southeastern Tibet. Geology. 2005;33:525–8.

    Google Scholar 

  86. Wang E, Kirby E, Furlong KP, van Soest M, Xu G, Shi X, Kamp PJJ, Hodges KV. Two-phase growth of high topography in eastern Tibet during the Cenozoic. Nat Geosci. 2012;5:640–5.

    Google Scholar 

  87. Wang P, Scherler D, Liu-Zeng J, Mey J, Avouac JP, Zhang Y, Shi D. Tectonic control of Yarlung Tsangpo gorge revealed by a buried canyon in southern Tibet. Science. 2014;346:978–81.

    CAS  PubMed  Google Scholar 

  88. Xing Y, Ree RH. Uplift-driven diversification in the Hengduan Mountains, a temperate biodiversity hotspot. Proc Natl Acad Sci. 2017;114:E3444–51.

    CAS  PubMed  PubMed Central  Google Scholar 

  89. Guo P, Zhao E. Pareas stanleyi—a record new to Sichuan, China and a key to the Chinese species. Asiatic Herpetological Research. 2004;10:280–1.

    Google Scholar 

  90. Yan F, Zhou W, Zhao H, Yuan Z, Wang Y, Jiang K, Jin J, Murphy RW, Che J, Zhang Y. Geological events play a larger role than Pleistocene climatic fluctuations in driving the genetic structure of Quasipaa boulengeri (Anura: Dicroglossidae). Mol Ecol. 2013;22:1120–33.

    PubMed  Google Scholar 

  91. Sun Y, Hu H, Huang H, Vargas-Mendoza CF. Chloroplast diversity and population differentiation of Castanopsis fargesii (Fagaceae): a dominant tree species in evergreen broad-leaved forest of subtropical China. Tree Genet Genomes. 2014;10:1531–9.

    Google Scholar 

  92. Mittermeier RA, Turner WR, Larsen FW, Brooks TM, Gascon C. Global biodiversity conservation: the critical role of hotspots. In: Zachos FE, Habel JC, editors. Biodiversity hotspots: distribution and protection of priority conservation areas. Berlin: Springer; 2011. p. 3–22.

    Google Scholar 

  93. Igea J, Tanentzap AJ. Multiple macroevolutionary routes to becoming a biodiversity hotspot. Science Advances. 2019;5:eaau8067.

    CAS  PubMed  PubMed Central  Google Scholar 

  94. Xiong JL, Gu HJ, Gong TJ, Zeng XM. Redescription of an enigmatic salamander, Pseudohynobius puxiongensis (Fei et Ye, 2000) (Urodela: Hynobiidae). Zootaxa. 2011;2919:51–9.

    Google Scholar 

  95. Yuan S, Barrett SCH, Duan T, Qian X, Shi M, Zhang D. Ecological correlates and genetic consequences of evolutionary transitions from distyly to homostyly. Ann Bot. 2017;120:775–89.

    PubMed  PubMed Central  Google Scholar 

  96. Wang W, Wei S, Chen M, Wu H. Female-biased dispersal of the Emei moustache toad (Leptobrachium boringii) under local resource competition. Asian Herpetological Research. 2019;10:24–31.

    Google Scholar 

  97. Bono JM, Pigage HK, Wettstein PJ, Prosser SA, Pigage JC. Genome-wide markers reveal a complex evolutionary history involving divergence and introgression in the Abert’s squirrel (Sciurus aberti) species group. BMC Evol Biol. 2018;18:139.

    PubMed  PubMed Central  Google Scholar 

  98. Hinojosa JC, Koubínová D, Szenteczki MA, Pitteloud C, Dincă V, Alvarez N, Vila R. A mirage of cryptic species: genomics uncover striking mitonuclear discordance in the butterfly Thymelicus sylvestris. Mol Ecol. 2019;28:3857–68.

    PubMed  Google Scholar 

  99. Lambert SM, Streicher JW, Fisher-Reid MC, Méndez de la Cruz FR, Martínez-Méndez N, García-Vázquez UO, Montes de Oca AN, Wiens JJ. Inferring introgression using RADseq and DFOIL: power and pitfalls revealed in a case study of spiny lizards (Sceloporus). Mol Ecol Resour. 2019;19:818–37.

    PubMed  Google Scholar 

  100. Larkin MA, Blackshields G, Brown NP, Chenna R, McGettigan PA, McWilliam H, et al. Clustal W and Clustal X version 2.0. Bioinformatics. 2007;23:2947–8.

  101. Rozas J, Ferrer-Mata A, Sánchez-DelBarrio JC, Guirao-Rico S, Librado P, Ramos-Onsins SE, Sánchez-Gracia A. DnaSP 6: DNA sequence polymorphism analysis of large data sets. Mol Biol Evol. 2017;34:3299–302.

    CAS  PubMed  Google Scholar 

  102. Bouckaert R, Heled J, Kühnert D, Vaughan T, Wu CH, Xie D, Suchard MA, Rambaut A, Drummond AJ. BEAST 2: a software platform for Bayesian evolutionary analysis. PLoS Comput Biol. 2014;10:e1003537.

    PubMed  PubMed Central  Google Scholar 

  103. Lanfear R, Frandsen PB, Wright AM, Senfeld T, Calcott B. PartitionFinder 2: new methods for selecting partitioned models of evolution for molecular and morphological phylogenetic analyses. Mol Biol Evol. 2017;34:772–3.

    CAS  PubMed  Google Scholar 

  104. Macey JR, Schulte JA, Larson A, Fang Z, Wang Y, Tuniyev BS, Papenfuss TJ. Phylogenetic relationships of toads in the Bufo bufo species group from the eastern escarpment of the Tibetan Plateau: a case of vicariance and dispersal. Mol Phylogenet Evol. 1998;9:80–7.

  105. Hoffman EA, Blouin MS. Evolutionary history of the northern leopard frog: reconstruction of phylogeny, phylogeography, and historical changes in population demography from mitochondrial DNA. Evolution. 2004;58:145–59.

    PubMed  Google Scholar 

  106. Sanguila MB, Siler CD, Diesmos AC, Nuneza O, Brown RM. Phylogeography, geographic structure, genetic variation, and potential species boundaries in Philippine slender toads. Mol Phylogenet Evol. 2011;61:333–50.

    PubMed  Google Scholar 

  107. Rambaut A, Suchard MA, Xie W, Drummond AJ. Tracer v1.6. 2014. Available from http://beast.bio.ed.ac.uk/Tracer.

  108. Sun X, Liu D, Zhang X, Li W, Liu H, Hong W, et al. SLAF-seq: an efficient method of large-scale de novo SNP discovery and genotyping using high-throughput sequencing. PLoS One. 2013;8:e58700.

    CAS  PubMed  PubMed Central  Google Scholar 

  109. Andrews S. FastQC: a quality control tool for high throughput sequence data. 2010. Retrieved from http://www.bioinformatics.babraham.ac.uk/projects/fastqc/.

  110. Catchen J, Hohenlohe PA, Bassham S, Amores A, Cresko WA. Stacks: an analysis tool set for population genomics. Mol Ecol. 2013;22:3124–40.

    PubMed  PubMed Central  Google Scholar 

  111. Paris JR, Stevens JR, Catchen JM. Lost in parameter space: a road map for stacks. Methods Ecol Evol. 2017;8:1360–73.

    Google Scholar 

  112. Danecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA, et al. The variant call format and VCFtools. Bioinformatics. 2011;27:2156–8.

    CAS  PubMed  PubMed Central  Google Scholar 

  113. Excoffier L, Lischer HEL. Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Mol Ecol Resour. 2010;10:564–7.

  114. Chang CC, Chow CC, Tellier LC, Vattikuti S, Purcell SM, Lee JJ. Second-generation PLINK: rising to the challenge of larger and richer datasets. GigaScience. 2015;4:7.

    PubMed  PubMed Central  Google Scholar 

  115. Peakall R, Smouse PE. GenAlEx 6.5: genetic analysis in Excel. Population genetic software for teaching and research-an update. Bioinformatics. 2012;28:2537–9.

  116. Felsenstein J. PHYLIP (phylogeny inference package) version 3.695. Department of Genome Sciences and Department of biology, University of Washington. 2013. Retrieved from http://evolution.genetics.washington.edu/phylip/.

  117. Pickrell JK, Pritchard JK. Inference of population splits and mixtures from genome-wide allele frequency data. PLoS Genet. 2012;8:e1002967.

    CAS  PubMed  PubMed Central  Google Scholar 

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

    Google Scholar 

  119. Phillips SJ, Anderson RP, Dudík M, Schapire RE, Blair ME. Opening the black box: an open-source release of Maxent. Ecography. 2017;40:887–93.

    Google Scholar 

  120. Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A. Very high resolution interpolated climate surfaces for global land areas. Int J Climatol. 2005;25:1965–78.

    Google Scholar 

  121. Synes NW, Osborne PE. Choice of predictor variables as a source of uncertainty in continental-scale species distribution modelling under climate change. Glob Ecol Biogeogr. 2011;20:904–14.

    Google Scholar 

  122. Dormann CF, Elith J, Bacher S, Buchmann C, Carl G, Carré G, et al. Collinearity: a review of methods to deal with it and a simulation study evaluating their performance. Ecography. 2013;36:27–46.

    Google Scholar 

  123. Araújo MB, New M. Ensemble forecasting of species distributions. Trends Ecol Evol. 2007;22:42–7.

    PubMed  Google Scholar 

  124. Phillips SJ, Dudík M. Modeling of species distributions with Maxent: new extensions and a comprehensive evaluation. Ecography. 2008;31:161–75.

    Google Scholar 

  125. Ledo RMD, Colli GR. The historical connections between the Amazon and the Atlantic Forest revisited. J Biogeogr. 2017;44:2551–63.

    Google Scholar 

  126. Gent PR, Danabasoglu G, Donner LJ, Holland MM, Hunke EC, Jayne SR, et al. The community climate system model version 4. J Clim. 2011;24:4973–91.

    Google Scholar 

  127. Watanabe S, Hajima T, Sudo K, Nagashima T, Takemura T, Okajima H, et al. MIROC-ESM 2010: model description and basic results of CMIP5-20c3m experiments. Geosci Model Dev. 2011;4:845–72.

    Google Scholar 

Download references

Acknowledgements

We would like to thank Guannan Wen and Yundong Gao for their helpful advice and discussions.

Funding

The design of the study, data collection and analysis, interpretation of data, and manuscript preparation were supported by the Ministry of Education of China through the Foundation of Key Laboratory of Southwest China Wildlife Resources Conservation (XNYB16–3), the Ministry of Science and Technology of China through the National Key Programme of Research and Development Project (2017YFC0505202), and the National Natural Science Foundation of China (31372181 and 31572243).

Author information

Authors and Affiliations

  1. Chengdu Institute of Biology, Chinese Academy of Sciences, #9 of Section 4, Ren-Min-Nan Road, Wuhou District, Chengdu, 610041, Sichuan Province, China

    Yuchi Zheng, Qiang Dai, Xianguang Guo & Xiaomao Zeng

  2. Key Laboratory of Southwest China Wildlife Resources Conservation (Ministry of Education), China West Normal University, Nanchong, 637009, Sichuan Province, China

    Yuchi Zheng

Authors
  1. Yuchi Zheng
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Qiang Dai
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Xianguang Guo
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Xiaomao Zeng
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

YZ and XZ conceived the project and conducted field sampling. YZ obtained the molecular data. YZ and QD analyzed the data. YZ, QD, XG, and XZ designed and wrote the manuscript. All authors have read and approved the manuscript.

Corresponding author

Correspondence to Yuchi Zheng.

Ethics declarations

Ethics approval and consent to participate

This work was carried out with the permission of the Ethics Committee of the Chengdu Institute of Biology (20181064C).

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

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

Supplementary information

Additional file 1: Table S1.

Sample information and haplotype designation. (XLS 87 kb)

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. 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 in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Zheng, Y., Dai, Q., Guo, X. et al. Dynamics behind disjunct distribution, hotspot-edge refugia, and discordant RADseq/mtDNA variability: insights from the Emei mustache toad. BMC Evol Biol 20, 111 (2020). https://0-doi-org.brum.beds.ac.uk/10.1186/s12862-020-01675-8

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://0-doi-org.brum.beds.ac.uk/10.1186/s12862-020-01675-8

Keywords

BMC Ecology and Evolution

ISSN: 2730-7182

Contact us