Mitochondrial characteristics of Pseudohynobius flavomaculatus a protected salamander in China, and biogeographical implications for the family Hynobiidae (Amphibia, Caudata)

Pseudohynobius flavomaculatus a provincially-protected salamander species, inhabits mountainous areas of Chongqing and sur-rounding 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 orig inated 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.


Introduction
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 (http://www.iucnredlist.org/). 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.

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 (8 th 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 HiSeq TM 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 (https://github.com/ndierckx/NOVOPlasty; (Dierckxsens et al. 2016).

Gene annotation and sequence analysis
Assembled mitochondrial genes of P. flavomaculatus were annotated using the Mitos web server (http://mi-tos2.bioinf.uni-leipzig.de/index.py; (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 (http://lowelab.ucsc.edu/ tRNAscan-SE/). 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).
We looked up each taxon's risk of extinction based on the IUCN Red List of Threatened Species (http://www. iucnredlist.org/). 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 (http://tree.bio.ed.ac. uk/software/figtree/).

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).
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 . 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).

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).

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.  (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). 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 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. 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. 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. 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 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 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 (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.

Conclusion
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.