Research Article
Research Article
Mitochondrial characteristics of Pseudohynobius flavomaculatus a protected salamander in China, and biogeographical implications for the family Hynobiidae (Amphibia, Caudata)
expand article infoYu Zhang, Meng Wang, Ruli Cheng, Yang Luo, Yingwen Li, Zhihao Liu, Qiliang Chen, Yanjun Shen
‡ Chongqing Normal University, Chongqing, China
Open Access


Pseudohynobius flavomaculatus a provincially-protected salamander species, inhabits mountainous areas of Chongqing and surrounding provinces in China. In the present study, the complete mitochondrial genome of P. flavomaculatus was sequenced and analyzed. The mitogenome is 16,401 bp in length and consisted of 13 protein-coding genes, 2 ribosomal RNA genes, 22 transfer RNA genes, and a control region. We performed a novel phylogenetic analysis, which demonstrated a sister relationship between P. flavomaculatus and P. jinfo. The 95% confidence interval around our new divergence date estimate suggest that Hynobiidae originated at 101.62–119.84 (mean=110.87) Ma. Species within Hynobiidae diverged successively in the Cenozoic era, and hynobiid speciation coincides primarily with geologic events. Our biogeographical inference demonstrates that nearly all early hynobiids divergences correspond to geological estimates of orogeny, which may have contributed to the notably high dN/dS ratio in this clade. We conclude that orogeny is likely a primary, dynamic factor, which may have repeatedly initiated the process of speciation in the family Hynobiidae.

Key Words

biogeography, divergence time, mitochondrial genome, molecular phylogenetics


Salamanders of the family Hynobiidae are an abundant component of the amphibian fauna endemic to Asia, and exhibit particularly high density of species in East Asia, where they inhabit montane distributions from 100 to 4,400 m (Yang et al. 2008; Xiong et al. 2010). The clade represents one of the first-diverging groups of extant terrestrial vertebrates; thus, studies of the family are of particular interest to researchers exploring the origin and evolution of terrestrial vertebrates (Zhao and Hu 1984). The discovery of many new hynobiid species has greatly enriched our understanding of diversification in the family. However, several controversies exist regarding taxonomic description of species, such as Batrachuperus tibetanus, Liua tsinpaensis, and Pseudohynobius flavomaculatus (Xiong et al. 2009). Among these unresolved questions, the validity of the genus Pseudohynobius, and the taxonomic validity of P. flavomaculatus, has been disputed since the first reports of these taxa (Fei and Ye 1982; Fei and Ye 1983); some named taxa were not convincingly validated until more than two decades later (Fei and Ye 1982; Fei and Ye 1983; Zhao and Wu 1995; Xu 2002; Li et al. 2004; Zeng et al. 2006).

Subsequent studies involving the family Hynobiidae family revealed apparently-resolved phylogenetic relationships (Zhang et al. 2006; Li et al. 2011; Zheng et al. 2011; Tang et al. 2015). Zhang et al.’s (2006) study suggested that phylogenetic relationships among living species could best be understood in the context of their geographical distribution; these authors interpreted hynobiids as divisible into five groups, based on the region each inhabits, namely Central Asia, East, and North Asia, and West versus, Central China (Zhang et al. 2006). Whereas the geographic distribution of hynobiids has been interpreted as mainly influenced by Cenozoic tectonic evolution (Zhang et al. 2006; Li et al. 2011), studies have shown that the hynobiids originated in the Middle Cretaceous (Zhang et al. 2006; Li et al. 2011; Zheng et al. 2011; Tang et al. 2015). To date, only one plausible hypothesis exists for the center of origin, of the clade: “out of North China” (Zhang et al. 2006). This interpretation is contingent on a 161 Ma-year-old fossil Chunerpeton tianyiensis (originating in Ningcheng County, Inner Mongolia), which is the earliest-known salamander fossil (Gao and Shubin 2003; Zhang et al. 2006).

Pseudohynobius flavomaculatus is distributed in Sichuan, Chongqing, Guizhou, Hunan, and Hubei Provinces of China. It inhabits moss or soil cavities under roots of bamboo and shrubs, in high-altitude mountains (Yang et al. 2008). Excessive hunting and overall habitat degradation by human activities have led to a noticeable reduction in hynobiid populations (Zhou et al. 2004; Yang et al. 2008). Thus, in 2019, P. flavomaculatus was grouped as a vulnerable (VU) species in the IUCN Red List of Threatened Species ( To ensure protection of P. flavomaculatus, it was listed as a key, provincial protected animal in the Chongqing municipality of China.

The mitochondrion is a vital, semiautonomous organelle, which is present in most eukaryotic cells. It exhibits maternal inheritance and is inherited through the cytoplasm of the ovum (Avise and Saunders 1984; Xia and Wang 1998; Boore 1999). Mitochondrial DNA of multicellular animals consists of a small, closed, circular molecule, with a rapid evolutionary rate of molecular evolution; it does not undergo recombination and has a size range of approximately 15–20 kilobases (kb; (Boore 1999; Dai et al. 2013). The composition of the mitogenome is relatively constant, comprising 13 protein-coding genes (PCGs), two ribosomal RNA (rRNA) genes, 22 transfer RNA (tRNA) genes, and a control region (CR; (Wolstenholme 1992; Boore 1999). Therefore, because it has been so well characterized, mitochondrial DNA has been extensively employed for studying evolutionary genetics of animals, genetic structure and diversity of populations, molecular ecology, and species identification (Shen et al. 2017; Besenbacher et al. 2019; Su and Liang 2019; Shen et al. 2020; Kvist et al. 2022; Zuo et al. 2022).

In this study, we assembled the complete mitogenome of P. flavomaculatus and analyzed the mitogenomic structure, nucleotide composition, codon usage, and tRNA secondary structure. Furthermore, we estimate phylogenetic relationships, characterize rates of synonymous versus nonsynonymous substitutions (dN/dS ratios), and infer a divergence timeframe for Hynobiidae, based on 13 PCGs. Finally, we provide a quantitative biogeographical analysis of rage evolution in Hynobiidae, based on novel findings presented in this study.

Materials and methods

Sample collection and DNA extraction

In this study, the P. flavomaculatus specimens were obtained from Pengshui County (29°15'N, 108°11'E), Chongqing municipality of China, and were humanely sacrificed in the field, while adhering to the recommended laws and regulations. The specimens were stored in anhydrous ethanol in the Chongqing Key Laboratory of Animal Biology, Chongqing Normal University. Total DNA was extracted from the dehydrated muscle tissues by using the TaKaRa MiniBEST Universal Genomic DNA Extraction Kit Version 5.0 (TaKaRa Biotech, Beijing, China). Extracted DNA was stored at –20 °C. All animal experiments were performed according to the Guideline for the Care and Use of Laboratory Animals (8th Edn., National Academies Press) and approved by the Committee of Laboratory Animal Experimentation at Chongqing Normal University (Approval No. Zhao-2021-713-02).

Mitogenome sequencing and assembly

The extracted mitogenomic DNA was visualized via 1% agarose gel electrophoresis, and the target band was excised and gel-extracted. DNA was physically fragmented, into a range of 300–500-bp, using ultrasound, and was then purified to construct a sequence library. Library construction involved the following steps: DNA terminal repair, addition of an adaptor to the 3’ end, adaptor ligation, recovery of the target fragment through agarose gel electrophoresis, and PCR amplification of the final target fragment. The mitogenome was sequenced using Illumina HiSeqTM and raw sequence reads were processed using the following steps: (1) removal of the adapter sequence from reads, (2) removal of bases containing non–AGCT at the 5’ end, (3) trimming off the tail end of low-quality reads (Phred quality score less than Q20), (4) trimming off of the reads containing 10% ambiguous N bases, and (5) disposal of sequences less than 50 bp in length. Clean data were de novo assembled using NOVOPlasty software (; (Dierckxsens et al. 2016).

Gene annotation and sequence analysis

Assembled mitochondrial genes of P. flavomaculatus were annotated using the Mitos web server (http://mitos2.bioinf.uni–; (Bernt et al. 2013). The graphical map of the mitogenome was drawn using the visualization tool Circos (Krzywinski et al. 2009) and modified manually. The secondary structures of tRNA genes were estimated using the Mitos web server and the tRNA-Scan SE web server ( Codon usage was determined using MEGA 7.0 (Kumar et al. 2016). Start and stop codons of the PCGs were determined, based on the vertebrate mitochondrial genetic code table (Osawa et al. 1989). Strand asymmetry was calculated using the following formulae: AT-skew = (A – T)/(A + T); GC-skew = (G – C)/(G + C) (Perna and Kocher 1995).

Phylogenetic analyses

To reconstruct the phylogenetic relationships of species of the family Hynobiidae, we downloaded the complete mitogenomes of 31 hynobiids and four outgroup species from GenBank (; see Suppl. material 3: Table S1). These 31 hynobiid species included two species of Pseudohynobius, six Batrachuperus, 14 Hynobius, two Liua, three Onychodactylus, two Salamandrella, one Ranodon, and one Pachyhynobius. Four outgroups included two species of Cryptobranchidae (Andrias davidianus, A. japonicus) and two species of Salamandridae (Echinotriton andersoni, E. chinhaiensis).

We looked up each taxon’s risk of extinction based on the IUCN Red List of Threatened Species ( MEGA 7.0 (default parameter settings) was used to align each PCG of 36 mitogenomes. Finally, 13 PCGs, with 11,085-nt length, were arranged 5’ to 3’ .

The concatenated mitogenome (sequences of 13 PCGs) was used for phylogenetic analyses, which were performed using Bayesian inference (BI) and the Maximum Likelihood (ML) topology estimation methods. The best-fit model GTR + G + I was selected for both ML and BI analyses through jModeltest (Darriba et al. 2012). Our Maximum Likelihood analysis was conducted using PhyML 3.0, with 1,000 bootstrap replicates, under out best-fit model (Mccormack et al. 2009), whereas the BI analyses were performed using MrBayes 3.2.7, with 5,000,000 generations (Ronquist and Huelsenbeck 2003), and sampled trees every 1,000 generations. We conservatively discarded as burn-in (Altekar et al. 2004) the initial 50% cycles from our Markov Chain Monte Carlo (MCMC) iterations. Convergence was assessed by evaluating split frequencies, and was considered to be achieved when the average standard deviation of split frequencies reached <0.01.

Divergence time estimation

The temporal framework for hynobiid divergence was estimated using Bayesian analyses in BEAST version 2.6.0, and based on 13 concatenated PCGs. A BEAST XML file was generated in BEAUTi 2, using a model consisting of the GTR, with the gamma parameter set to 4. Two calibration methods were used in this study: (1) the divergence time between B. yenyuanensis and other Batrachuperus species (~22.6–26 Ma; (Zhang et al. 2006), and (2) Onychodactylus diverging from its closest relatives at 106.4–114.9 Ma; (Zhang et al. 2006; Tang et al. 2015). Our MCMC chains of 50,000,000 generations were sampled every 1,000 generations, and we discarded a conservative initial 50% samples as burn-in. The ESSs in TRACER version 1.7.1 (Rambaut et al. 2018) were used to evaluate the convergence of each parameter (ESSs > 200). The sample tree was annotated using TreeAnnotator and then visualized (not using a fixed tree topology) in FigTree version 1.4.3 (

dN/dS ratios

The dN/dS ratios of 31 Hynobiidae species were analyzed using EasyCodeML version 1.21 (Gao et al. 2019) based on 13 concatenated PCGs. Ratios of dN/dS were based on the preferred ML phylogenetic tree (Fig. 3) and branch lengths were estimated using a free-ratio (model = 1) branch model. The branch lengths were estimated using a free-ratio (model = 1) branch model. We set parameters to κ=3, ω=1.3, and α=0, respectively.

Results and discussion

Mitogenome organization and nucleotide composition

The complete mitogenome of P. flavomaculatus is 16,401 bp in length and comprises 22 tRNAs (one for each amino acid, two for leucine, and two for serine), two rRNA genes (rrnL and rrnS), 13 PCGs (cox1-3, nad1-6, nad4l, cob, atp6, and atp8), and one CR (Fig. 1, Table 1) (GenBank accession: MT476485). Of the 37 genes, 28 were located on the heavy (+) strand, whereas the remaining nine genes, namely trn (Q, A, N, C, Y, S2, E, P) and nad6, were located on the light (–) strand. The size overlap between nine gene pairs was 39 nt, with the longest overlap of 15 bp between nad5 and nad6. The mitogenome has nine intergenic regions of 175 bp size, with the longest region of 128 bp between trnT and trnP, and the size of remaining intergenic regions ranging from 1 to 33 nt (Table 1). The base composition of the P. flavomaculatus mitogenome was: A = 33.45%, T = 31.83%, C = 20.60%, and G = 14.11%, indicating bias towards AT nucleotides (65.28%); atp8 was the most distinct of all the genes, with the AT content of 72.62% (Table 2). However, compared with the AT nucleotide percentage of other hynobiids (range: 63.3%–68.8%), that of P. flavomaculatus is not conspicuously anomalous (Suppl. material 3: Table S2). Of the 13 PCGs, both the AT-skew (only cox2, atp8, and nad2 were positive) and GC-skew (only nad6 was positive) of total PCGs were negative (–0.0481 and –0.1884, respectively), indicating that T and C have a higher proportion than A and G (Table 2). The AT-skew of the complete mitogenome of P. flavomaculatus was positive (0.0248), whereas the GC-skew (only nad6 and tRNAs were positive) was negative (–0.1869). We also determined the base composition of complete mitogenomes of all other species of Hynobiidae; their AT-skew was slightly positive, ranging 0.00 to 0.04, whereas the GC-skew was negative, ranging from –0.21 to –0.17 (Suppl. material 3: Table S2).

Table 1.

Summary of the P. flavomaculatus mitogenome; TA and T are incomplete termination codons.

Gene Direction Location Size Anticodon Amino acid (length) Start codon Stop codon Intergenic nucletide
trnF + 1–69 69 TTC
rrnS + 70–1,002 933 –1
trnV + 1,002–1,071 70 GTA
rrnL + 1,072–2,673 1,602
trnL2 + 2,674–2,748 75 TTA
nad1 + 2,749–3,719 971 323 ATG TA
trnI + 3,720–3,790 71 ATC –1
trnQ 3,790–3,861 72 CAA –1
trnM + 3,861–3,929 69 ATG
nad2 + 3,930–4,973 1,044 347 ATG TAA –2
trnW + 4,972–5,040 69 TGA 2
trnA 5,043–5,111 69 GCA 1
trnN 5,113–5,185 73 AAC 33
trnC 5,219–5,285 67 TGC
trnY 5,286–5,353 68 TAC 1
cox1 + 5,355–6,905 1,551 516 GTG TAA
trnS2 6,906–6,976 71 TCA 2
trnD + 6,979–7,047 69 GAC 2
cox2 + 7,050–7,737 688 229 ATG T
trnK + 7,738–7,808 71 AAA 1
atp8 + 7,810–7,977 168 55 ATG TAA –10
atp6 + 7,968–8,651 684 227 ATG TAA –1
cox3 + 8,651–9,435 785 261 ATG TA –1
trnG + 9,435–9,502 68 GGA
nad3 + 9,503–9,851 349 116 GTG T
trnR + 9,852–9,920 69 CGA
nad4l + 9,921–10,217 297 98 ATG TAA –7
nad4 + 10,211–11,588 1,378 459 ATG T
trnH + 11,589–11,657 69 CAC
trnS1 + 11,658–11,725 68 AGC
trnL1 + 11,726–11,797 72 CTA
nad5 + 11,798–13,618 1,821 606 ATG TAA –15
nad6 13,604–14,122 519 172 ATG AGA
trnE 14,123–14,191 69 GAA 5
cob + 14,197–15,337 1,141 380 ATG T
trnT + 15,338–15,406 69 ACA 128
trnP 15,535–15,606 72 CCA
CR + 15,607–16,401 795
Table 2.

Composition and skewness of the P. flavomaculatus mitogenome.

Region Size(bp) A(%) T(%) G(%) C(%) A+T(%) AT-skew GC-skew
Mitogenome 16,401 33.45 31.83 14.11 20.60 65.28 0.03 –0.19
nad1 971 31.51 33.37 13.08 22.04 64.88 –0.03 –0.26
nad2 1,044 34.58 34.20 10.54 20.69 68.78 0.01 –0.33
cox1 1,551 27.98 34.11 17.67 20.25 62.09 –0.10 –0.07
cox2 688 33.14 32.56 14.53 19.77 65.70 0.01 –0.15
atp8 168 38.69 33.93 10.12 17.26 72.62 0.07 –0.26
atp6 684 32.46 33.63 11.99 21.93 66.09 –0.02 –0.29
cox3 785 29.55 33.50 15.92 21.02 63.05 –0.06 –0.14
nad3 349 31.23 34.38 13.18 21.20 65.61 –0.05 –0.23
nad4l 297 30.98 36.70 11.45 20.88 67.68 –0.09 –0.29
nad4 1,378 33.60 33.82 11.90 20.68 67.42 –0.00 –0.27
nad5 1,821 33.00 34.49 12.36 20.15 67.49 –0.02 –0.24
nad6 519 19.65 43.74 25.43 11.18 63.39 –0.38 0.39
cob 1,141 29.27 32.78 14.20 23.66 62.05 –0.06 –0.25
PCGs 11,396 31.14 34.29 14.02 20.53 65.43 –0.05 –0.19
tRNAs 1,539 33.72 31.97 18.26 16.05 65.69 0.03 0.06
rRNAs 2,535 38.97 26.27 16.37 18.38 65.24 0.20 –0.06
CR 795 29.90 35.00 16.70 18.40 64.90 –0.08 –0.05
Figure 1. 

Gene map of the P. flavomaculatus mitogenome. The internal gray circulars are GC content graph.

PCGs and codon usage

The 13 PCGs of P. flavomaculatus contained 7 NADH genes (nad1-6 and nad4l), 2 ATP genes (atp6 and atp8), and four cytochrome genes (cox1-3 and cob) (Fig. 1, Table 1). Among the 13 PCGs, only one NADH gene (nad6) was located on the light (–) strand, whereas the remaining 12 PCGs were located on the heavy (+) strand. The start codon of 11 PCGs was ATG, whereas that of two PCGs (cox1, nad3) was GTG. The stop codon for the gene nad6 was AGA, that for six PCGs (nad2, nad4l, nad5, cox1, atp8, atp6) was TAA, that for two PCGs (nad1 and cox3) was TA, and that for four PCGs (cox2, nad3, nad4, cob) was single T. The incomplete stop codons (TA and T) were presumably completed as TAA following post-transcriptional polyadenylation (Anderson et al. 1981). In addition, termination of the gene nad6 by AGA and other characteristics are similar to those of the PCGs of other species of hynobiid salamanders (Zhang et al. 2003; Zhou et al. 2004; Lee et al. 2011; Tang et al. 2015; Huang et al. 2016).

Table 3 and Fig. 2 show the relative synonymous codon usage (RSCU) of the 13 PCGs of P. flavomaculatus mitogenome. The 13 PCGs consisted of 3,695 codons, with no start or stop codon. Among the amino acids deduced from the codons, Leu (16.32%) had the highest translation frequency, whereas Cys (0.78%) had the lowest translation frequency. In addition, the third codon position with the maximum RSCU value of each amino acid codon type was found to be generally occupied by A or U, for example, Ile (AUU), Phe (UUU), Arg (CGA), and Val (GUA) (Table 3). This finding suggests that A and T are used more frequently in the nucleotide sequences of PCGs than C and G. This phenomenon is consistent with the A + T bias (65.43%) in the nucleotide composition of 13 PCGs. In most metazoan mitogenomes, synonymous substitutions usually occur at the third codon position (Zhang et al. 2013). The higher representation of nucleotides A and T at the third codon position may be associated with genome bias, optimal selection of tRNA, or the DNA repair efficiency (Crozier and Crozier 1993; Chai and Du 2012; Fischer et al. 2013). Moreover, the higher representation of the nucleotides A and T leads to a subsequent bias in the encoded amino acids (Shen et al. 2012). In the final analysis, this pattern of codon usage may be affected by mutation pressure and natural selection (Uddin and Chakraborty 2019).

Table 3.

The codon number and relative synonymous codon usage in P. flavomaculatus mitochondrial protein-coding genes.

Codon Count RSCU Codon Count RSCU Codon Count RSCU Codon Count RSCU
UUU(F) 199 1.66 UCU(S) 45 0.97 UAU(Y) 82 1.43 UGU(C) 15 1.03
UUC(F) 41 0.34 UCC(S) 47 1.01 UAC(Y) 33 0.57 UGC(C) 14 0.97
UUA(L) 267 2.66 UCA(S) 132 2.85 UAA(*) 0 0.00 UGA(W) 91 1.70
UUG(L) 29 0.29 UCG(S) 8 0.17 UAG(*) 0 0.00 UGG(W) 16 0.30
CUU(L) 110 1.09 CCU(P) 37 0.77 CAU(H) 61 1.27 CGU(R) 17 0.97
CUC(L) 32 0.32 CCC(P) 27 0.57 CAC(H) 35 0.73 CGC(R) 7 0.40
CUA(L) 143 1.42 CCA(P) 112 2.35 CAA(Q) 82 1.80 CGA(R) 41 2.34
CUG(L) 22 0.22 CCG(P) 15 0.31 CAG(Q) 9 0.20 CGG(R) 5 0.29
AUU(I) 249 1.55 ACU(T) 76 1.16 AAU(N) 85 1.12 AGU(S) 22 0.47
AUC(I) 72 0.45 ACC(T) 54 0.83 AAC(N) 67 0.88 AGC(S) 24 0.52
AUA(M) 194 1.67 ACA(T) 122 1.87 AAA(K) 78 1.77 AGA(*) 0 0.00
AUG(M) 38 0.33 ACG(T) 9 0.14 AAG(K) 10 0.23 AGG(*) 0 0.00
GUU(V) 53 1.16 GCU(A) 65 1.01 GAU(D) 43 1.30 GGU(G) 37 0.68
GUC(V) 22 0.48 GCC(A) 58 0.90 GAC(D) 23 0.70 GGC(G) 41 0.75
GUA(V) 88 1.93 GCA(A) 117 1.81 GAA(E) 74 1.54 GGA(G) 91 1.66
GUG(V) 19 0.42 GCG(A) 18 0.28 GAG(E) 22 0.46 GGG(G) 50 0.91
Figure 2. 

Relative synonymous codon usage in the P. flavomaculatus mitogenome.

tRNAs and rRNAs

The 22 tRNA genes ranged in length from 67 to 75 nt. Of these, 14 genes were located on the heavy (+) strand, and eight genes were located on the light (–) strand (Fig. 1, Table 1). The estimation results of tRNA secondary structures by the two web servers were consistent. The analysis of secondary structures of 22 tRNAs showed that all the tRNAs had a typical cloverleaf structure, except trnS1, which was t-shaped (Suppl. material 1: Fig. S1).

The two rRNA genes had a total length of 2,535 nt (AT% = 65.24, AT-skew = 0.1947), with the rrnS (12S rRNA) gene being 933-nt long and the rrnL (16S rRNA) gene being 1,602-nt long (Table 2). Both rRNA genes were located on the heavy (+) strand. The rrnL and rrnS genes were separated by trnV and located between trnL2 and trnF (Fig. 1, Table 1).

Phylogenetic analysis

As depicted in Fig. 3, the results of the two trees from BI and ML analyses were mostly consistent, and both divided the species of the family Hynobiidae into five clades. The slight topological structure differences between BI and ML trees reflect two primary differences. The first topological structure difference resulted from the phylogenetic location of H. yiwuensis. In Clade I, Hynobius yiwuensis and H. amjiensis grouped together in our ML analysis (ML bootstrap: BP = 65); however, H. yiwuensis and a clade involving H. chinensis, H. guabangshanensis, and H. maoershanensis clustered together in our BI analysis (BI posterior probabilities: PP = 0.53). Although a few studies have investigated the phylogenetic relationship of H. yiwuensis, most earlier studies (Li et al. 2011; Zheng et al. 2011; Tang et al. 2015) are in agreement with results of our BI analysis. The second major difference involves Clade II (Fig. 3) which contains Sublades II-1, II-2, II-3, and II-4. Subclade II-2 and Clade II-1 clustered together in our ML analysis (BP = 66), however Clade II-2 and Clade II-3 clustered were most closely-related in our BI analysis (PP = 1). Previous studies have reported that Clade II-2 (Liua) is more closely related to Clade II-3 (Pseudohynobius) than to Clade II-1 (Batrachuperus) (Peng et al. 2010; Zheng et al. 2011; Tang et al. 2015).

Figure 3. 

Phylogeny of P. flavomaculatus based on nucleotide sequences of 13 protein-coding genes determined using the BI and ML methods. Numbers next to the nodes indicate bootstrap (ML, site left of “/”) and posterior probability (BI, right). The two separate branches on the right belong to the BI tree, which indicate two differences in topology with the ML tree.

Both ML and BI phylogenetic trees indicated that the genus Hynobius is a monophyletic group (ML bootstrap: BP = 100; BI posterior probabilities: PP = 1; Fig. 3). The genera Pseudohynobius, Liua, and Salamandrella and some species of Batrachuperus were grouped into Clade II. The genus Pachyhynobius was sister to all hynobiids (Claded I & II). The genus Ranodon and two species of the genus Batrachuperus were grouped into Clade IV. The genus Onychodactylus was sister to all other hynobiids and is most likely the first-diverging lineage. In Clade II-3, P. flavomaculatus and P. jinfo clustered in a single node with high nodal support values (ML bootstrap: BP = 100; BI posterior probabilities: PP = 1; Fig. 3), indicating that P. flavomaculatus and P. jinfo exhibit a sister relationship.

The phylogenetic tree indicated six species of the genus Batrachuperus, of which four were clustered into Clade II-1, and the remaining two species (B. mustersi and B. gorganensis) and Ranodon sibiricus were clustered into Clade IV. These results are consistent with phylogenetic relationships reconstructed in previous analyses, although B. mustersi and B. gorganensis have been classified under the genus Paradactylodon in earlier studies (Peng et al. 2010; Tang et al. 2015). Batrachuperus mustersi and B. gorganensis have been referred to as Paradactylodon in some other studies as well, and the nomenclature of the genus Paradactylodon has been debated for many years (Stock et al. 2019); in one earlier study, Paradactylodon was considered a junior synonym for Batrachuperus (Reilly 1987). Another study highlighted that Batrachuperus is diphyletic, comprising a Central–Western Asia group and a West China group (Thorn and Raffaëlli 2001). Therefore, the convergent (or symplesiomorphic) phenotypic characteristics of these two groups may be attributed to similar selective environments or retention of ancestral characters (Zhang et al. 2006). Despite the morphological support for this genus, the Central–Western Asia group does not exhibit a close phylogenetic relationship with the West China group. Conversely, the Central–Western Asia group, comprising B. mustersi and B. gorganensis, appears more closely-related to Ranodon (Fig. 3).

Divergence time estimation

Divergence date estimates of within Hynobiidae are shown in Fig. 4. The molecular dating results showed that the family-level clade appeared in the Jurassic (162.79 Ma, 95% confidence interval: 145.41–186.37 Ma). The geological date estimated for the fossil Chunerpeton tianyiense (~161 Ma) indicates that the divergence between Cryptobranchidae and Hynobiidae occurred in Asia before the Middle Jurassic, with a minimally conservative division at approximately 145 Ma (Gao and Shubin 2003). An earlier study showed that the diversification of Hynobiidae began at approximately 110.7 Ma (confidence interval: 106.4–114.9 Ma) (Zhang et al. 2006), which is fairly consistent with the results of (Tang et al. 2015). In this study, we estimated 110.87 Ma (confidence interval: 101.62–119.84 Ma), which is additionally consistent with both previous studies.

Figure 4. 

Divergence time estimation of the Hynobiidae. Node bars indicate 95% CIs of the divergence time estimate. Chronostratigraphic scale data are derived from International Commission on Stratigraphy (

The initial hynobiid divergences, approximately 110–112 Ma, gave rise to the first-branching genus (Onychodactylus) of the family Hynobiidae. Over the next 40 million years, the existing branches of the family evolved gradually. At the Paleogene (Cenozoic era), results of our analyses suggests that salamanders of the Hynobiidae began to diversify again. Divergence time estimates suggest that the genus Hynobius, split from other hynobiids of Clade II at 59.31 Ma (confidence interval: 52.09–66.46 Ma). Coincidentally, the ancestors of the Central–Western Asia group and West China group of Batrachuperus appeared at approximately the same time period (i.e., 45 Ma; Fig. 4).

Hynobiids are widely distributed over Asia, with Clades I and V distributed mainly in East Asia, Clades II and III in Western and Central China, and Clade IV in Central Asia. Hynobiids have low vagility, and their dispersal is limited by mountains and deserts (Zhang et al. 2006) . According to Zhang et al.’s (2006) hypothesis, living hynobiids originated in northern China in the Middle Cretaceous when the continent was relatively flat and wet, and geological conditions facilitated dispersal of hynobiids over broad geographical ranges. Since 50–60 Ma, mountains such as the Himalayas and the Kunlun ranges have been forced upwards as a result of the collision between the Eurasian plate and the Indian plate (Fang 2017). We assume that mountains act as barriers to inland dispersal of hynobiids, thereby creating geographical isolation, and contributing to diversification in allopatry (Zhang et al. 2006). Consequently, the lack of genetic communication among hynobiid lineages, during lengthy periods of isolation, may have led to the emergence of different species (Shen et al. 2020). Furthermore, speciation and geographical distribution of hynobiids in coastal and other regions appear closely related to Cenozoic tectonic evolution (Zhang et al. 2006; Li et al. 2011; Tang et al. 2015).

dN/dS ratios and biogeography analysis

The dN/dS ratio (ω) is the ratio of nonsynonymous (dN) to synonymous (dS) substitutions rates (ω = dN/dS). Values of ω > 1, ω = 0, and ω = < 1 indicate that the effect of natural selection on the change of PCGs is positive, neutral, and negative, respectively (Kryazhimskiy and Plotkin 2008). In the present study, the dN/dS ratios of PCGs for all species were <1 (Fig. 5). The maximum dN/dS ratio of 0.4242 was noted for Pachyhynobius shangchengensis, whereas the minimum dN/dS ratio of 0.0337 was noted for H. guabangshanensis. In addition, these two most closely related species exhibited similar dN/dS ratios.

Figure 5. 

The dN/dS ratios of the 13 concatenated PCGs for 32 Hynobiidae species. The Roman number represents the corresponding partition in the phylogenetic tree.

Interestingly, one finding was surprising and conceptually intriguing: We found a systematic relationship between timing of diversification and dN/dS ratios: the earlier the divergence (Fig. 4), the higher the dN/dS ratios (Fig. 5). Because of the notable rates of diversification in hynobiids occurring in the Cenozoic (Fig. 4) and frequent geological activity of mountains in Central and East Asia during this era (Wu 2003; Wang et al. 2014; Zhang 2015; Deng 2016), we find it reasonable to speculate that geological processes may have influenced rates of synonymous and nonsynonymous substitution in hynobiids.

Among the living hynobiids, P. shangchengensis appeared at the earliest time (approximately 60 Ma) and was distributed in the Dabie Mountains (located at the junction of Anhui, Hubei, and Henan Provinces) (Zhao et al. 2013). The Dabie Mountains had the highest uplift scale between 70 and 60 Ma (Wu 2003). Thus, we speculate that the P. shangchengensis speciated is likely associated with the pronounced, frequent environmental changes, which characterize these mountain ranges (Edelaar 2018). Most of the hynobiids in Clade II appeared approximately 33.5–15 Ma, except for P. flavomaculatus and P. jinfo, which diverged 4.01 Ma. Both species are distributed in the Wuling Mountains, and the geological activity that created this massif occurred before the Cretaceous period, with no significant changes in the Cenozoic (Xiao 2017). In Clade II, the lineage leading to P. puxiongensis diverged earliest (~33 Ma), whereas B. yenyuanensis appeared slightly later (~24 Ma); both these species are distributed in the Daliang Mountains areaSichuan Province. The Daliang Mountains, located on the southeast edge of the Tibetan Plateau, was rapidly uplifted during 40–20 Ma and rapidly created a landform, similar to the present day (Deng 2016). In addition, B. pinchonii (15.69 Ma) and B. tibetanus (17.35 Ma) are widely distributed in the east of Tibet Plateau, and B. londongensis (15.69 Ma) is distributed in Mount Emei. The areas of Himalaya, Gangdise, and western Kunlun were uplifted rapidly during the period 20–13 Ma, and the Emei mountain experienced the fastest uplift during the period 20–10 Ma (Zhang et al. 2010; Zhou et al. 2020).

The genus Liua comprises only two species, namely L. tsinpaensis and L. shihi, which are distributed in the Qinling mountains and Daba Mountains, respectively. Liua tsinpaensis (22.27 Ma) was mainly distributed on the northern edge of the Qinling Mountain, which was gradually uplifted during 25–10 Ma and accompanied by frequent taphrogeny (Meng 2017). The Daba Mountains, where L. shihi inhabits, have been continuously uplifting since 153 Ma, and they experienced gradual uplift during 50–5 Ma (Zhang 2015).

Divergence time estimation indicated that H. chinensis and H. guabangshanensis (0.17 Ma), H. maoershanensis (5.92 Ma), H. arisanensis and H. formosanus (1.83 Ma), P. flavomaculatus and P. jinfo (4.01 Ma), and O. zhaoermii and O. fischeri (0.39 Ma) diverged later, and their dN/dS ratios are low. However, H. chinensis, H. guabangshanensis, H. maoershanensis, P. flavomaculatus, P. jinfo do not appear to have speciated in a manner associated with geologic events. Hynobius arisanensis & H. formosanus may have arisen as a result of geological folds in the central mountain range of Taiwan (Ye 1983), and O. zhaoermii & O. fischeri evolved possibly due to volcanic activity in the Changbai Mountain area (Liu et al. 2015).

Our data indicate that the earlier-diverging hynobiids, with high dN/dS ratios, may have experienced orogeny initially associated with their speciation, whereas the later hynobiids, with low dN/dS ratios, did not experience similar orogenic. This finding suggests that continuous environmental variation associated with topographical complexity may spur environmental selection on low-vagility terrestrial vertebrates like hynobiid salamanders, which appears associated with increased nonsynonymous substitutions of gene sequences (Shen et al. 2020). In contrast, hynobiids that had evolved in relatively stable environments retained fewer nonsynonymous substitutions during the same period. Therefore, the dN/dS ratios of earlier hynobiids were high, whereas those of later hynobiids were low.


In this study, we sequenced and annotated the complete mitogenome of P. flavomaculatus. This mitogenome was 16,401 nt in length and contained 37 genes and one CR. The start codon of 11 PCGs was ATG, whereas that of two PCGs was GTG; the stop codons for 6, 2, 4, and 1 PCGs were TAA, TA, T, and AGA, respectively. The nucleotide composition of the complete mitogenome was biased towards AT nucleotides. Divergence time estimation indicated that the Hynobiidae originated at 110.87 Ma. Additionally, the dN/dS ratios of earlier-diverging lineages of hynobiids were found to be higher than those of later-diverging hynobiids. Together with our biogeographic interpretation, we suggest that earlier hynobiids may have higher dN/dS ratios because of frequent and dramatic environmental changes, whereas the later-diverging lineages may have low dN/dS ratios because of a relatively stable environmental set of conditions through time.

Conflicts of interest

The authors declare that they have no conflict of interest.


The authors sincerely thank all the crews for their help with the manuscript writing and data analysis. This work was supported by the Science and Technology Research Program of Chongqing Municipal Education Commission (Grant No. KJQN202100503).


  • Altekar G, Dwarkadas S, Huelsenbeck JP, Ronquist F (2004) Parallel metropolis coupled Markov chain Monte Carlo for Bayesian phylogenetic inference. Bioinformatics 20(3): 407–415.
  • Anderson S, Bankier AT, Barrell BG, de Bruijn MHL, Coulson AR, Drouin J, Eperon IC, Nierlich DP, Roe BA, Sanger F, Schreier PH, Smith AJH, Staden R, Young IG (1981) Sequence and organization of the human mitochondrial genome. Nature 290(5806): 457–465.
  • Avise JC, Saunders NC (1984) Hybridization and introgression among species of sundish (Lepomis): Analysis by mitochondrial DNA and allozyme markers. Genetics 108(1): 237–255.
  • Bernt M, Donath A, Juhling F, Externbrink F, Florentz C, Fritzsch G, Pütz J, Middendorf M, Stadler PF (2013) MITOS: Improved de novo metazoan mitochondrial genome annotation. Molecular Phylogenetics and Evolution 69(2): 313–319.
  • Besenbacher S, Hvilsom C, Marques-Bonet T, Mailund T, Schierup MH (2019) Direct estimation of mutations in great apes reconciles phylogenetic dating. Nature Ecology & Evolution 3(2): 286–292.
  • Chai H, Du Y (2012) The complete mitochondrial genome of the pink stem borer, Sesamia inferens, in comparison with four other Noctuid moths. International Journal of Molecular Sciences 13(8): 10236–10256.
  • Dai L, Zhu B, Liu Q, Wei G, Liu C (2013) Characterization of the complete mitochondrial genome of Bombyx mori strain H9 (Lepidoptera: Bombycidae). Gene 519(2): 326–334.
  • Darriba D, Taboada GL, Doallo R, Posada D (2012) jModelTest 2: More models, new heuristics and parallel computing. Nature Methods 9(8): 772–772.
  • Deng B (2016) Cenozoic mountain-building processes in the Daliangshan, southeastern margin of the Tibetan Plateau: Evidence from low-temperature thermochronology and thermal modeling. Chinese Journal of Geophysics-chinese Edition 59(6): 2162–2175
  • Dierckxsens N, Mardulyn P, Smits G (2016) NOVOPlasty: De novo assembly of organelle genomes from whole genome data. Nucleic Acids Research 45(4): e18.
  • Fang X (2017) Phased uplift of the Tibetan Plateau. Science & Technology Review 35(6): 42–50.
  • Fei L, Ye CY (1982) A new species of salamander from Hubei province–Hynobius flavomaculatus. Dong Wu Fen Lei Xue Bao 7(2): 225–228.
  • Fei L, Ye CY (1983) Systematic discussion of the genus Hynobiidae with description of a new species. Acta Zootax Sin 2(4): 31–37.
  • Fischer C, Koblmuller S, Gully C, Schlotterer C, Sturmbauer C, Thallinger GG (2013) Complete mitochondrial DNA sequences of the threadfin cichlid (Petrochromis trewavasae) and the blunthead cichlid (Tropheus moorii) and patterns of mitochondrial genome evolution in cichlid fishes. PLoS ONE 8(6): e67048.
  • Gao F, Chen C, Arab DA, Du Z, He Y, Ho SYW (2019) EasyCodeML: A visual tool for analysis of selection using CodeML. Ecology and Evolution 9(7): 3891–3898.
  • Krzywinski M, Schein J, Birol I, Connors J, Gascoyne R, Horsman D, Jones SJ, Marra MA (2009) Circos: An information aesthetic for comparative genomics. Genome Research 19(9): 1639–1645.
  • Kumar S, Stecher G, Tamura K (2016) MEGA7: Molecular Evolutionary Genetics Analysis version 7.0 for bigger datasets. Molecular Biology and Evolution 33(7): 1870–1874.
  • Kvist S, Utevsky S, Marrone F, Ben Ahmed R, Gajda Ł, Grosser C, Huseynov M, Jueg U, Khomenko A, Oceguera-Figueroa A, Pešić V, Pupins M, Rouag R, Sağlam N, Świątek P, Trontelj P, Vecchioni L, Müller C (2022) Extensive sampling sheds light on species-level diversity in Palearctic Placobdella (Annelida: Clitellata: Glossiphoniiformes). Hydrobiologia 849(5): 1–21.
  • Lee B, Kim JY, Song S, Hur JM, Cho JY, Park YC (2011) The complete mitochondrial genome sequence of the Kori salamander Hynobius yangi (Caudata: Hynobiidae). Mitochondrial DNA 22(5–6): 168–170.
  • Li Y, Wu M, Wang XL (2004) Phylogenetic relationships of Hynobiidae based on sequences of mitochondrial 16S ribosomal RNA gene. Dong Wu Xue Bao 50(3): 464–469.
  • Liu JQ, Chen SS, Guo WF, Sun CQ, Zhang ML, Guo ZF (2015) Research advances in the Mt. Changbai volcano. Bulletin of Mineralogy. Petrology and Geochemistry 34(4): 710–723.
  • Mccormack JE, Huang H, Knowles LL (2009) Maximum Likelihood Estimates of Species Trees: How Accuracy of Phylogenetic Inference Depends upon the Divergence History and Sampling Design. Systematic Biology 58(5): 501–508.
  • Osawa S, Ohama T, Jukes TH, Watanabe K (1989) Evolution of the mitochondrial genetic code I. Origin of AGR serine and stop codons in metazoan mitochondria. Journal of Molecular Evolution 29(3): 202–207.
  • Peng R, Zhang P, Xiong J, Gu H, Zeng X, Zou F (2010) Rediscovery of Protohynobius puxiongensis (Caudata: Hynobiidae) and its phylogenetic position based on complete mitochondrial genomes. Molecular Phylogenetics and Evolution 56(1): 252–258.
  • Perna NT, Kocher T (1995) Patterns of nucleotide composition at fourfold degenerate sites of animal mitochondrial genomes. Journal of Molecular Evolution 41(3): 353–358.
  • Rambaut A, Drummond AJ, Xie D, Baele G, Marc A (2018) Posterior Summarization in Bayesian Phylogenetics Using Tracer 1.7. Systematic Biology 67(5): 901–904.
  • Shen X, Li X, Sha Z, Yan B, Xu Q (2012) Complete mitochondrial genome of the Japanese snapping shrimp Alpheus japonicus (Crustacea: Decapoda: Caridea): Gene rearrangement and phylogeny within Caridea. Science China-life Sciences 55(7): 591–598.
  • Shen Y, Kou Q, Zhong Z, Li X, He L, He S, Gan X (2017) The first complete mitogenome of the South China deep‐sea giant isopod Bathynomus sp. (Crustacea: Isopoda: Cirolanidae) allows insights into the early mitogenomic evolution of isopods. Ecology and Evolution 7(6): 1869–1881.
  • Stock M, Fakharzadeh F, Kuhl H, Rozenblutkościsty B, Leinweber S, Patel RP, Forster DW (2019) Shedding Light on a Secretive Tertiary Urodelean Relict: Hynobiid Salamanders (Paradactylodon persicus s.l.) from Iran, Illuminated by Phylogeographic, Developmental, and Transcriptomic Data. Genes 10(4): 306.
  • Su T, Liang A (2019) Comparative analysis of seven mitochondrial genomes of Phymatostetha (Hemiptera: Cercopidae) and phylogenetic implications. International Journal of Biological Macromolecules 125: 1112–1117.
  • Thorn R, Raffaëlli J (2001) Les Salamandres de l’Ancien Monde. Paris: Boubée, 380 pp.
  • Wu Q (2003) The uplift process research of post-orgeny action in Dabie area. Chinese Academy of Sciences Guangzhou.
  • Xia D, Wang W (1998) The study of animal mitochondrial DNA and its application to the study on fish population genetic structure. Shuichan Xuebao 22(4): 364–370.
  • Xiao H (2017) The geological study on shale gas reservoirs of Longmaxi formation in east Sichuan Wuling mountain area. China University of Geosciences.
  • Xiong JL, Sun P, Zhu WW, Liu XY (2009) Review on the Taxonomic Variation of Genus and Species in Chinese Hynobiids. Sichuan Journal of Zoology 28(6): 958–961.
  • Xiong JL, Zhu WW, Sun P (2010) The systematic taxonomy of Chinese endemic genera in Hynobiidae. Agricultural Science and Technology 11(1): 55–58.
  • Xu J (2002) Some taxonomic problems of the Hynobiidae. Acta Scientiarum Naturalium Universitatis Sunyatseni 41(1): 79–81.
  • Yang L, Gong D, Mu M (2008) Present research situation and resource conservation of Hynobiidae (Amphibia: Urodela) in China. Shengtaixue Zazhi 27(1): 111–116.
  • Ye S (1983) BiologiaTectonic characteristics of Taiwan and its adjacent sea areas. Tropic Oceanology 2(4): 260–268.
  • Zeng X, Fu J, Chen L, Tian Y, Chen X (2006) Cryptic species and systematics of the hynobiid salamanders of the LiuaPseudohynobius complex: Molecular and phylogenetic perspectives. Biochemical Systematics and Ecology 34(6): 467–477.
  • Zhang Y (2015) Study of tectonic thermal evolution history and uplift in Mesozoic-cenozoic in the north rim of Sichuan basin. Chang’an University.
  • Zhang P, Chen Y, Zhou H, Wang X, Qu L (2003) The complete mitochondrial genome of a relic salamander, Ranodon sibiricus (Amphibia: Caudata) and implications for amphibian phylogeny. Molecular Phylogenetics and Evolution 28(3): 620–626.
  • Zhang P, Chen Y, Zhou H, Liu Y, Wang X, Papenfuss TJ, Wake DB, Qu L-H (2006) Phylogeny, evolution, and biogeography of Asiatic Salamanders (Hynobiidae). Proceedings of the National Academy of Sciences of the United States of America 103(19): 7360–7365.
  • Zhang K, Wang G, Ji J, Luo M, Kou X, Wang Y, Xu YD, Chen FN, Chen RM, Song BW, Zhang JY, Liang YP (2010) Paleogene-Neogene stratigraphic realm and sedimentary sequence of the Qinghai-Tibet Plateau and their response to uplift of the plateau. Science China-earth Sciences 53(9): 1271–1294.
  • Zhao EM, Hu QX (1984) Studies on Chinese Tailed Amphibians. Sichuan Science and Technology Press.
  • Zhao EM, Wu GF (1995) Taxonomic status of Ranodon tsinpaensis Liu & Hu, 1966, with a discussion on Pseudohynobius flavomaculatus (Fei & Ye, 1982) as its synonym. Sichuan Journal of Zoology 14(1): 20–23.
  • Zhao Y, Zhang Y, Li X (2013) Molecular Phylogeography and Population Genetic Structure of an Endangered Species Pachyhynobius shangchengensis (hynobiid Salamander) in a Fragmented Habitat of Southeastern China. PLoS ONE 8(10): e78064.
  • Zheng Y, Peng R, Kuro-o M, Zeng X (2011) Exploring Patterns and Extent of Bias in Estimating Divergence Time from Mitochondrial DNA Sequence Data in a Particular Clade: A Case Study of Salamanders (Order Caudata). Molecular Biology and Evolution 28(9): 2521–2535.
  • Zhou Z, Wu J, Wang GZ, Deng JH, Deng B, Luo Q et al. (2020) Cenozoic accelerated erosion of the Emeishan, eastern margin of Tibetan Plateau. Chinese Journal of Geophysics-chinese Edition 63(4): 1370–1385.
  • Zuo Q, Zhang Z, Shen Y (2022) Novel mitochondrial gene rearrangements pattern in the millipede Polydesmus sp. GZCS‐2019 and phylogenetic analysis of the Myriapoda. Ecology and Evolution 12(3): e8764.

Supplementary materials

Supplementary material 1 

Figure S1

Yu Zhang, Meng Wang, Ruli Cheng, Yang Luo, Yingwen Li, Zhihao Liu, Qiliang Chen, Yanjun Shen

Data type: JPG file

Explanation note: Secondary structures of the 22 transfer RNA (tRNA) genes of P. flavomaculatus. tRNAs are labelled with the abbreviations of their corresponding amino acids. The lines between the English letters indicate the phosphodiester bonds and Watson–Crick pairing.

This dataset is made available under the Open Database License ( The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited.
Download file (741.85 kb)
Supplementary material 2 

Figure S2

Yu Zhang, Meng Wang, Ruli Cheng, Yang Luo, Yingwen Li, Zhihao Liu, Qiliang Chen, Yanjun Shen

Data type: JPG file

Explanation note: The dN/dS ratios of the each PCG for 32 Hynobiidae species.

This dataset is made available under the Open Database License ( The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited.
Download file (329.65 kb)
Supplementary material 3 

Tables S1, S2

Yu Zhang, Meng Wang, Ruli Cheng, Yang Luo, Yingwen Li, Zhihao Liu, Qiliang Chen, Yanjun Shen

Data type: Adobe PDF file

Explanation note: Table S1. List of Caudata species analyzed in this study with their GenBank accession numbers, and endangered category from the IUCN Red List of Threatened Species. Table S2. Composition and skewness of the mitogenome for 36 Caudata species.

This dataset is made available under the Open Database License ( The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited.
Download file (417.74 kb)