A single juvenile specimen of the genus Pseudohynobius (Fig. 2A, B) was collected in August 2024 from Xishui County, Guizhou Province, China (Fig. 1). The habitat was a pond in a bamboo forest with a large number of subadults; no eggs or adults were found. Over the following near two years, 16 field surveys were conducted, but no additional adults or juveniles were observed. In addition, we also detected the presence of a Pseudohynobius population in the Fanjingshan Mountains, Guizhou Province, China, although no specimen was obtained (Fig. 2C). All experimental procedures strictly adhered to the National Standard Guidelines of the People’s Republic of China (GB/T 35892-2018). Euthanasia was humanely performed using tricaine methanesulfonate (MS-222), and the specimen was subsequently preserved in the Animal Ecology Laboratory at Guizhou Normal University, Guiyang, Guizhou Province, China.

Figure 2. 

Two cryptic species of the genus Pseudohynobius. A, B. Specimen from Xishui County, Guizhou, China; C. Specimen photographed in Fanjingshan National Nature Reserve, Guizhou, China, by Mr. Xian-Wei Meng, who authorized its use.

To clarify the phylogenetic placement of the juvenile Pseudohynobius specimen collected from Xishui County, Guizhou Province, China, we sequenced its mitochondrial genome. Total genomic DNA was extracted from each sample from 95% ethanol-preserved tissues using the cetyltrimethylammonium bromide method. Illumina sequencing libraries with 300–500 bp insert fragments were generated for each sample and sequenced on the Illumina NovaSeq 6000 sequencer at TSING KE Biological Technology Co. (Chengdu, China), generating approximately 12 Gb of raw data. The raw data were filtered using fastp v.0.23.4 (Chen 2023) using default parameters, after which the mitogenome was assembled using MitoZ v.3.6 (Meng et al. 2019).

In this study, we used a total of 24 mitochondrial sequences from the genus Pseudohynobius for phylogenetic analysis, including a newly sequenced mitochondrial genome, as well as five mitochondrial genomes, 15 cytochrome b (Cyt b) sequences, two COI sequences, and two 16S rRNA sequences downloaded from GenBank (Table 1) (Fu et al. 2001; Zeng et al. 2006; Zhang et al. 2006; Peng et al. 2010; Xia et al. 2012; Yang et al. 2013; Huang et al. 2016; Zhang et al. 2022). We constructed two sequence matrices: dataset 1, which contains only Cyt b and is used for species delimitation and calculation of genetic distances; and dataset 2, which includes complete mitogenomes and concatenated sequences of some species (Cyt b, COI, and 16S rRNA) for assessing divergence times and inferring biogeography. The mitogenome was extracted using PhyloSuite v.1.2.3 (Zhang et al. 2020a) for 13 protein-coding genes, two rRNAs, and 22 tRNAs for subsequent genetic analysis.

Multiple sequence alignment was performed using MAFFT v.7.4 (Katoh and Standley 2013) within PhyloSuite v.1.2.3 (Zhang et al. 2020a) and checked using MEGA v.7.0 (Kumar et al. 2016) to rule out possible errors. In addition, PartitionFinder v.2.1.1 (Lanfear et al. 2017) was used to select the best-fit partitioning and nucleotide substitution model for the two datasets based on the Bayesian information criterion. In dataset 1, each gene fragment was preset as an independent partition. We followed a previous study and selected the genus Liua (two species), the genus Batrachuperus (four species), and Salamandrella keyserlingii as outgroups (Table 1).

For dataset 1, phylogenetic trees were reconstructed using Bayesian inference (BI) and maximum likelihood (ML) methods based on best-fit partitioning and nucleotide substitution models. Bayesian inference analysis was performed using MrBayes v.3.2.1 (Ronquist et al. 2012). Each BI analysis was run independently using four Markov chain Monte Carlo chains (three heated chains and one cold chain) starting with a random tree; each chain was run for two million generations and sampled every 1,000 generations. Convergence of the data runs was confirmed when the average standard deviation of split frequencies was less than 0.01. The ML analysis was performed using IQ-TREE v.2.0.4 (Nguyen et al. 2015) based on the best-fit model with 10,000 ultrafast bootstrap (UFB) replicates (Hoang et al. 2018). The ML analysis was performed until a correlation coefficient of at least 0.99 was achieved.

To infer whether the juvenile sample represents a cryptic species and to identify possible cryptic species within the genus Pseudohynobius, we performed species delimitation. We employed three species delimitation methods based on the concept of phylogenetic species: Assemble Species by Automatic Partitioning (ASAP) (Puillandre et al. 2021), based on sequence variation; the Bayesian implementation of the Poisson Tree Processes (bPTP) (Zhang et al. 2013) and Multi-rate Poisson Tree Processes (mPTP) (Kapli et al. 2017), both based on phylogenetic trees. Tree-based approaches incorporate phylogenetic relationships and are particularly useful when branch lengths reflect evolutionary divergence, while ASAP focuses on pairwise genetic distances and tends to be more conservative (Zhang et al. 2013; Kapli et al. 2017; Puillandre et al. 2021). This complementary use of multiple methods allows for a more comprehensive and reliable assessment, especially in light of potential variation in sampling strategy and population history. For ASAP (https://bioinfo.mnhn.fr/abi/public/asap, accessed January 11, 2025), the calculations were performed online using dataset 1 and the default parameters. For bPTP (https://species.h-its.org, accessed January 11, 2025) and mPTP (http://mptp.h-its.org, accessed January 11, 2025), the calculations were performed online using the ML tree and the time tree of dataset 1, with the default parameters. Final species definition results were obtained from the consensus of the three delimitation methods.

Genetic distances were calculated using the uncorrected p-distance model with 1,000 bootstrap replications in MEGA v.7.0 (Kumar et al. 2016).

Molecular dating and time tree reconstruction were conducted in BEAST v.2.4.7 (Bouckaert et al. 2014), using dataset 2. We set two calibration nodes with reference to Chen et al. (2015): (1) the divergence of Salamandrella keyserlingii from Liua and Batrachuperus occurred at 33.3 Ma, with a 95% confidence interval (CI) of 38.6–28.4 Ma; (2) the MRCA of the genus Pseudohynobius was at 6.7 Ma (95% CI: 8.7–5.0 Ma). The BEAST analysis used an uncorrelated lognormal relaxed clock and a Yule tree prior. The analysis was run for 200 million generations and sampled every 1,000 generations. All calibrations used a normal prior, monophyly, and standard deviation values of 2.7 (for calibration node 1) and 1.0 (for calibration node 2). Convergence of run parameters was checked using Tracer v.1.7.1 (Rambaut et al. 2018) to ensure that the effective sample size of all parameters was greater than 200. Three runs were performed, and the resulting trees were merged using LogCombiner v.2.4.7 (Bouckaert et al. 2014). A maximum clade credibility tree was generated using TreeAnnotator v.2.4.7 (Bouckaert et al. 2014) with a burn-in of 25%.

The possible ancestral areas and dispersal routes of the genus Pseudohynobius were reconstructed using the package BioGeoBEARS (Matzke 2013). Prior to analysis, six models were evaluated to select the best-fitting model (Matzke 2014), including Dispersal–Extinction–Cladogenesis (DEC), a maximum likelihood version of Dispersal–Vicariance Analysis (DIVALIKE), and a Bayesian biogeographic inference model (BAYAREALIKE), each with and without founder-event speciation (+J). The maximum number of areas in the ancestral ranges was limited to two.

To infer the divergence dynamics of the genus Pseudohynobius, we followed the methods outlined in a previous phylogenetic study (Li et al. 2024). Results from BEAST were summarized to obtain time intervals with 95% confidence intervals for lineage divergence events. These events were defined as the maximum number of observed lineage divergence events (MDivE) per 0.5 Ma and were used to illustrate temporal trends. The data were smoothed by calculating means using a 0.5 Ma sliding window.

The total lengths of dataset 1 and dataset 2 were 1,140 bp (640 conserved sites, 500 variable sites, and 393 parsimony-informative sites) and 15,540 bp (10,147 conserved sites, 5,351 variable sites, and 3,509 parsimony-informative sites), respectively. For dataset 1, the nucleotide substitution models for the 1^st, 2^nd, and 3^rd codons are TVMEF+G, TRN+I, and HKY+G. The nucleotide substitution model for dataset 2 is shown in Table 2. The newly sequenced mitogenome was deposited in NCBI under accession number PV018879.

Species of the genus Pseudohynobius show a west-to-east distribution pattern in southwestern China, along the Hengduan Mountains, the Wumeng Mountains, the Dalou Mountains, and the Wuling–Wushan Mountains (Fig. 1). The overall trend in elevation decreases from west to east, with a mean elevation of about 1,750 m. Among them, P. puxiongensis, P. shuichengensis, and P. jinfo are distributed above 1,750 m, while P. flavomaculatus, P. guizhouensis, P. kuankuoshuiensis, and the two cryptic species occur below this elevation (Fig. 3).

Figure 3. 

Elevational ranges of species of the genus Pseudohynobius, with dashed lines indicating mean elevations.

Species delimitation based on 28 Cyt b sequences using three methods (ASAP, bPTP, and mPTP) assigned the sequences to eight phylogenetic species within Pseudohynobius (Suppl. materials 1–3). The newly collected single sample and P. flavomaculatus (DQ335719) were identified as independent phylogenetic species, designated Pseudohynobius sp. 1 and Pseudohynobius sp. 2, while P. flavomaculatus (MT476485) clustered with four samples of P. jinfo (Suppl. materials 1–3).

Both maximum likelihood and Bayesian analyses of dataset 1 produced the same topology (Fig. 4, Suppl. material 4). The BI and ML trees showed similar branch lengths and support values; within Pseudohynobius, the BI tree had 16 internal nodes with posterior probability support > 0.95 (Fig. 4A), and the ML tree had 12 internal nodes with bootstrap support > 95% (Fig. 4B). For species-level clades, low support was observed only between P. kuankuoshuiensis and two other species (Pseudohynobius sp. 1 and P. guizhouensis). Pseudohynobius can be divided into three clades: Clade I, the basal clade, includes only P. puxiongensis; Clade II contains P. flavomaculatus, Pseudohynobius sp. 2, and P. jinfo; and Clade III includes P. shuichengensis, P. kuankuoshuiensis, Pseudohynobius sp. 1, and P. guizhouensis (Fig. 4). A similar topology was obtained for the BI tree reconstructed from dataset 2 (Fig. 5), but it showed distinct inconsistencies between P. shuichengensis and P. kuankuoshuiensis (Suppl. material 5).

Figure 4. 

Phylogeny and species delimitation based on mitochondrial Cyt b. A. BI tree; Bayesian posterior probabilities (BPP) from BI analyses are shown at nodes; B. ML tree; ultrafast bootstrap supports (UFB) from ML analyses are shown at nodes. Branches are color-coded to represent different species based on the consensus results of species delimitation using ASAP, bPTP, and mPTP.

Genetic distances based on Cyt b revealed that the two cryptic lineages exhibit minimum interspecific distances of 3.2% and 4.2% compared to other congeners (Table 3), which are substantially higher than the recently reported minimum distance of 2.1% between Hynobius chinensis and H. guabangshanensis (Wang et al. 2023). This indicates that genetic differentiation between the two cryptic lineages and their closely related species has reached the species level.

Our molecular dating results suggest that Pseudohynobius originated in the middle Miocene, approximately 14.57 Ma (95% CI: 19.54–10.32) (node 1), and that the extant crown group originated in the late Miocene, approximately 8.62 Ma (95% CI: 10.48–6.80) (node 2) (Fig. 5A). Interspecific divergence occurred primarily between 8 Ma and 1.4 Ma (Fig. 5A).

Figure 5. 

Time tree, biogeographic history, and diversification dynamics based on mitogenomes. A. Time tree with ancestral area reconstruction. Blue bars at nodes indicate 95% CI of divergence times. Circled nodes represent the most probable ancestral regions; arrows show dispersal directions; B. Possible dispersal routes of the genus Pseudohynobius; C. Divergence dynamics based on the maximum number of observed lineage divergence events (MDivE) per Ma. OIHM: orogeny intensification of the Hengduan Mountains; D. Global temperature inferred from deep-sea oxygen isotope records (δ18O) in benthic foraminifera (Zachos et al. 2008), modeled mean annual precipitation under idealized CO₂ levels (Farnsworth et al. 2019).

Changes in clade divergence within Pseudohynobius were assessed over time using biogeographic dynamic meta-analysis. In total, nine clade divergence events were identified (Suppl. material 6). MDivE analyses showed that cladistic divergence began at approximately 10 Ma, increased rapidly between 8 and 5 Ma, peaked around 4 Ma, and declined sharply after 2 Ma (Fig. 5B).

The results of ancestral area reconstruction are shown in Fig. 5, Suppl. material 4. Models with “+J” had higher log-likelihood (LnL) values than those without, indicating a significantly better fit in all cases except for DIVALIKE (Table 4). Under the best-fit DEC+J model, vicariance and dispersal events likely shaped the current diversity and distribution of Pseudohynobius. Ancestral area reconstruction under this model suggests that the most recent common ancestor of extant Pseudohynobius most likely inhabited the Hengduan Mountains and the Dalou Mountains (Fig. 5A). An initial vicariant event separated the ancestral populations into Clade I and Clades II–III, after which the ancestral populations in the Dalou Mountains dispersed into the Wumeng Mountains, the Wuling–Wushan Mountains, and the Miaoling Mountains at 4.71 Ma, 3.19 Ma, and 1.44 Ma, respectively (Fig. 5A). In situ diversification in the Dalou and Wuling–Wushan Mountains may also have contributed to species formation, as in the cases of P. jinfo and P. flavomaculatus (Fig. 5A).

Our phylogeny is consistent with the results of Yang et al. (2013): P. puxiongensis represents the basalmost clade, Clade I; P. flavomaculatus and P. jinfo comprise Clade II; and P. guizhouensis, P. kuankuoshuiensis, and P. shuichengensis comprise Clade III. However, we also discovered that P. kuankuoshuiensis is only moderately supported as the sister clade of P. guizhouensis–Pseudohynobius sp. 1 (Figs 4, 5). DensiTree plots, in which the topologies of multiple time trees were superimposed, showed discordance between P. kuankuoshuiensis and P. shuichengensis (Suppl. material 5). We hypothesize that gene flow may be the main contributor to this discordance. Similar discordance was observed between Pseudohynobius, Liua, and Batrachuperus, which have a phylogenetic relationship of (Liua, (Batrachuperus, Pseudohynobius)) (Figs 4, 5) or (Pseudohynobius, (Liua, Batrachuperus)) (Zhang et al. 2022), but the nuclear gene tree supports a topology of (Batrachuperus, (Pseudohynobius, Liua)) (Chen et al. 2015). In addition, we also took into account the taxonomic controversy mentioned in the introduction, as it may influence the interpretation of phylogenetic and biogeographic reconstruction results. Evidence has been presented that gene flow, incomplete lineage sorting, and hybridization may all contribute to this discordance (Meleshko et al. 2021; Rivas-González et al. 2023), but the major contribution remains to be supported by phylogenomic evidence.

Integrating taxonomy for species delimitation is becoming a consensus, where phylogenetic species should be accepted when evidence of phylogenetic isolation is supported by other results, e.g., morphology, reproduction, and climatic niche (van Elst et al. 2025). In this study, phylogeny (Fig. 3), species delimitation (Suppl. materials 1, 2), and genetic distance (Table 3) support the existence of two cryptic phylogenetic species in Pseudohynobius: Pseudohynobius sp. 1 from Xishui, Guizhou (Fig. 2A, B), and Pseudohynobius sp. 2 from Tianping Mountain, Sangzhi County, Hunan. We also acknowledge that for the two identified cryptic lineages, species delimitation was based on Cyt b sequences from only a single sample, and no morphological comparisons were possible due to the lack of adult specimens and detailed morphological data. This limitation may introduce additional bias. Therefore, future efforts that expand field surveys to include more specimens, incorporate nuclear gene fragments or even genome-wide data, and obtain sufficient morphological evidence—particularly from adult individuals—will be essential to provide more robust support for species delimitation.

In addition, further surveys are needed in the Dalou Mountains, Wuling Mountains, and Wushan Mountains. For example, our 2022 photographs of unknown species of the genus Pseudohynobius in the Fanjingshan Mountains (Figs 1, 2C), combined with the two new cryptic species revealed in this study, suggest that the range and diversity of the genus Pseudohynobius need to be reassessed, which is valuable for formulating conservation policies and carrying out conservation actions.

For the mitogenome, the estimated origination time for the total group of Pseudohynobius (14.57 Ma vs. 49.85 Ma) and for crown group Pseudohynobius (8.62 Ma vs. 33.38 Ma) are younger than those reported by Zhang et al. (2022). The reasons for generating older divergence times are related to alternative saturation and the rapid evolution of mitochondrial DNA, resulting in rather ancient estimates (Zheng et al. 2011; Chen et al. 2015). However, we cannot simply assume that divergence times assessed using mitogenomes are always less precise than those derived from nuclear markers or that mitochondrial dates are always overestimated (Chen et al. 2015). When divergence times are less than 15 Ma, mitochondrial substitutional saturation is minimal, and estimates based on mitochondrial data can still be precise (Chen et al. 2015). Our origin and interspecific divergence times are in close agreement with the results of Chen et al. (2015) and are less than 15 Ma, so it can be assumed that it is reasonable for us to infer the timing of the origin and divergence of the genus Pseudohynobius.

Orogeny and low temperatures have promoted divergence and diversification within the genus Pseudohynobius. Orogeny creates a variety of environmental conditions, including climatic niches, new habitats, and dispersal barriers, while low temperatures may isolate populations in specific habitats—all of which are conducive to promoting speciation (Hoorn et al. 2013; Xing and Ree 2017; Rahbek et al. 2019a; Li et al. 2021). Our inferred MDivE curve (Fig. 5C) shows that the divergence rate of Pseudohynobius was high between 9 and 4 Ma, and the diversification trend coincides with the rapid uplift of the Hengduan Mountains and the Yunnan–Guizhou Plateau from the late Miocene to late Pliocene (Zhou and Chen 1993; Sun et al. 2011; Favre et al. 2015), suggesting that the diversification of Pseudohynobius may be closely linked to orogenic events in these regions.

As part of the continued expansion of orogeny—such as the uplift of the Qinghai–Tibet Plateau (Ding et al. 2022) and the right-lateral extrusion of Indochina (Zhang et al. 2009; Fyhn and Phach 2015)—the Hengduan Mountains and the Yunnan–Guizhou Plateau experienced rapid uplift during the late Miocene, peaking just before the late Pliocene (Zhou and Chen 1993; Zhang et al. 2009; Wang et al. 2012). Strong orogeny shaped high-elevation mountains and varied microhabitats, which may have driven differentiation of ecological niches among populations (Rahbek et al. 2019b). Topography and microclimatic changes may have led Pseudohynobius to favor narrow, high-altitude cold regions, thereby promoting adaptation and speciation at high elevations (Fig. 2) (Fei and Ye 2016). For instance, in P. jinfo, the capillaries are in close contact with and protrude toward the epidermis, while the dorsal pigmentation layer is more developed (Jing et al. 2017). Such features in the high-altitude frog Nanorana parkeri have been linked to adaptations to hypoxia and elevated ultraviolet radiation at high altitudes (Yang et al. 2019; Fu et al. 2022).

Topographic changes influence the formation and migration of rivers, which may also act as barriers to gene flow between distinct populations. At present, the Jinshajiang River, Wujiang River, and Qingshuijiang River serve as geographical boundaries separating the distributions of species within the genus (Fig. 1). U-Pb dating evidence based on detrital zircon grains suggests that the establishment of the modern Yangtze River occurred in the late Miocene, from about 11 to 5 Ma (Fu et al. 2021; He et al. 2023), which is consistent with the divergence of the basal clade P. puxiongensis from 10.48 to 6.8 Ma (Fig. 4). Interspecific divergence within Clade II, on the other hand, corresponds to the formation of the Wujiang and Miaoling Mountains during the Pliocene to early Pleistocene (Zhou and Chen 1993), which may have served as geographic barriers contributing to species divergence. Distribution patterns of batrachians inhabiting southwestern China are severely affected by orogeny and riverine barriers, with close examples including Batrachuperus yenyuanensis (Jia et al. 2019), Glyphoglossus yunnanensis (Zhang et al. 2020b), and the Paramesotriton caudopunctatus complex (Chen et al. 2011; Luo et al. 2021). In addition, the intensification of the Asian monsoon may have promoted plant growth and geomorphological erosion in southwestern China (Fig. 4C), thereby providing abundant food resources and geographic isolation (Antonelli et al. 2018; Ding et al. 2020), which in turn created new ecological and evolutionary opportunities for speciation.

Tao Luo, Jiang Zhou, and Huai-Qing Deng conceived and designed the research; Tao Luo, Zhong-Lian Wang, and Zi-Fa Zhao conducted field surveys and collected samples; Tao Luo and Ming-Yuan Xiao analyzed genetic data; and Tao Luo, Ning Xiao, and Jiang Zhou wrote and revised the manuscript. All authors have read and approved the final version of the manuscript.