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

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

Assembled mitochondrial genes of P. flavomaculatus were annotated using the Mitos web server (http://mitos2.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).

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 (https://www.ncbi.nlm.nih.gov/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 (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.

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

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.

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

Figure 1. 

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

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

Figure 2. 

Relative synonymous codon usage in the P. flavomaculatus mitogenome.

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

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 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 (https://stratigraphy.org/).

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

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.