Zoosyst. Evol. 98 (2) 2022, 263-274 | DOI 10.3897/zse.98.66578

> PENSUFT.

Ate
BERLIN

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

Yu Zhang’, Meng Wang!, Ruli Cheng!, Yang Luo’, Yingwen Li!, Zhihao Liu’,

Qiliang Chen!, Yanjun Shen?

1 Chongging Key Laboratory of Animal Biology, School of Life Sciences, Chongqing Normal University, Chongqing, 401331, China

https://zoobank. org/6013D51F-5A43-4E57-A750-BF31592CF287

Corresponding author: Yanjun Shen (shenyanjun@cqnu.edu.cn)

Academic editor: Rafe Brown Received 27 March 2021 Accepted 6 July 2022 Published 13 July 2022

Abstract

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.

Key Words

biogeography, divergence time, mitochondrial genome, molecular phylogenetics

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 etal. 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; Liet 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

Copyright Zhang, Y. et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use,
distribution, and reproduction in any medium, provided the original author and source are credited.

264 Zhang, Y. et al.: Mitogenome characteristics of P. flavomaculatus and biogeography of Hynobiidae

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 Sich-
uan, Chongqing, Guizhou, Hunan, and Hubei Provinces of
China. It inhabits moss or soil cavities under roots of bam-
boo 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 vul-
nerable (VU) species in the IUCN Red List of Threatened
Species (http://www.iucnredlist.org/). To ensure protec-
tion of P. flavomaculatus, it was listed as a key, provincial
protected animal in the Chongqing municipality of China.

The mitochondrion is a vital, semiautonomous or-
ganelle, which is present in most eukaryotic cells. It ex-
hibits 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, circu-
lar 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 mitoge-
nome is relatively constant, comprising 13 protein-cod-
ing 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 ge-
netics of animals, genetic structure and diversity of pop-
ulations, 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 struc-
ture, nucleotide composition, codon usage, and tRNA
secondary structure. Furthermore, we estimate phylo-
genetic 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.

zse.pensoft.net

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 1n anhy-
drous ethanol in the Chongqing Key Laboratory of Ani-
mal 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 Edn.,
National Academies Press) and approved by the Commit-
tee 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 ex-
cised 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 liga-
tion, recovery of the target fragment through agarose gel
electrophoresis, and PCR amplification of the final target
fragment. The mitogenome was sequenced using Illumi-
na HiSeq™ 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-qual-
ity reads (Phred quality score less than Q20), (4) trim-
ming 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 NOVOPIlas-
ty 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

Zoosyst. Evol. 98 (2) 2022, 263-274

the PCGs were determined, based on the vertebrate mito-
chondrial 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 (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 Cryptobrachidae (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 Max-
imum 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 conser-
vatively discarded as burn-in (Altekar et al. 2004) the
initial 50% cycles from our Markov Chain Monte Carlo
(MCMC) iterations. Convergence was assessed by evalu-
ating 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 oth-
er 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

265

al. 2015). Our MCMC chains of 50,000,000 generations
were sampled every 1,000 generations, and we discard-
ed 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 parame-
ter (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 K=3, o=1.3, and a=0, respectively.

Results and discussion

Mitogenome organization and nucleotide
composition

The complete mitogenome of P. flavomaculatus 1s
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 (cox/-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

zse.pensoft.net

266 Zhang, Y. et al.: Mitogenome characteristics of P. flavomaculatus and biogeography of Hynobiidae

qpeu

Pseudohynobuus flavomaculatus

Complete Mitochondrial Genome

Gicomplex I (NADH dehydrogenase)

complex III (ubichinol cytochrome c reductase)
Micomplex IV (cytochrome c oxidase)

ATP synthase

Mitransfer RNAs

WM ribosomal RNAs

trnF(ttc)

16,401 bp

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

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

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 (cox/-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.

zse.pensoft.net

The start codon of 11 PCGs was ATG, whereas that of
two PCGs (coxl, nad3) was GTG. The stop codon for
the gene nad6 was AGA, that for six PCGs (nad2, nad4l,
nad5, cox 1, atpS, 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 co-
don usage (RSCU) of the 13 PCGs of P. flavomaculatus
mitogenome. The 13 PCGs consisted of 3,695 codons,

Zoosyst. Evol. 98 (2) 2022, 263-274

267

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 REC:
rns + 70-1,002 933 -l
trnV + 1,002-1,071 70 GTA
rent + 1,072-2,673 1,602
trnL2 + 2,674-2,748 75 TTA
nad1 + 2,/49-3,719 971 323 ATG TA
trnl + 3,720-3,790 71 ATC -]
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 i
trnN - 5,113-5,185 73 AAC 33
trnC - 5,219-5,285 67 TGC
trnY - 5,286-5,353 68 TAC ih
cox + 5,355-6,905 155511 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 ai
trnK + 7,/38-7 808 71 AAA 1
atp8 + 7,810-7,977 168 5D ATG TAA -10
atp6 + 7,968-8,651 684 BRT 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
nad4| + 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 ae CTA
nad5 + 11,798-13,618 1,821 606 ATG TAA -15
nad6 - 13,604-14,122 519 INEZ 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

with no start or stop codon. Among the amino acids
deduced from the codons, Leu (16.32%) had the high-
est 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 compo-
sition of 13 PCGs. In most metazoan mitogenomes, syn-
onymous 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 selec-
tion (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 struc-
tures by the two web servers were consistent. The anal-
ysis 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
tRNA) gene being 933-nt long and the rvnL (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).

zse.pensoft.net

268 Zhang, Y. et al.: Mitogenome characteristics of P. flavomaculatus and biogeography of Hynobiidae

NOsa

Leu2 Ile

Gin Met Trp Ala Asn Cys Tyr

Ser2 Asp Lys Gly Arg’ His’ Serl Leul Glu Thr Pro

=_ uu = —<— — oo

Figure 2. Relative synonymous codon usage in the P. flavomaculatus mitogenome.

Table 2. Composition and skewness of the P. flavomaculatus
mitogenome.

Region Size(bp) A(%) T(%) G(%) C(%) A+T(%) AT- GC-
skew 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 = 331.23 34.38 13.18 21.20 65.61 -0.05 -0.23
nad4] 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

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

zse.pensoft.net

Clade II (Fig. 3) which contains Sublades II-1, H-2, H-3,
and II-4. Subclade II-2 and Clade II-1 clustered together
in our ML analysis (BP = 66), however Clade I-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 1s 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 hynobi-
ids (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 1s most likely the first-diverging lineage.
In Clade H-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 ex-
hibit 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

Zoosyst. Evol. 98 (2) 2022, 263-274

100
100

100

65

tol) 100/1

Hynobius chinensis
Hynobius guabangshanensis
Hynobius maoershanensis
Hynobius yiwuensis
Hynobius amijiensis
Aynobius unisacculus

100/1 Hynobius quelpaertensis
100/1 Aynobius yangi
Hynobius leechii

Hynobius nebulosus

Hynobius nigrescens
100/1 -Hynobius arisanensis
-Hynobius formosanus

-Aynobius kimur

66

Liua tsinpaensis
Liua shihi

100
84/0.78 Pseudohynobius jinfo

100/1 aA

-Batrachuperus musterst
‘Batrachuperus gorganensis
-Ranodon sibiricus

100/1

100/1

‘Andrias davidianus

100/1 Andrias japonicus
Echinotriton andersoni

Echinotriton chinhaiensis

Pseudohynobius flavomaculatus

Pseudohynobius puxiongensis
a)

Pachyhynobius shangchengensis

269

Aynobius maoershanensis
Hynobius yiwuensis

Hynobius chinensis
Hynobius guabangshanensis
-Hynobius amjiensis

L Batrachuperus pinchonii

Batrachuperus londongensis 1-1

——— Batrachuperus tibetanus

Batrachuperus yenyuanensis

Lina tsinpaensis

Liua shihi | 2 | l
Pseudohynobius flavomaculatus

Pseudohynobius jinfo

Pseudohynebius puxiongensis

Outgroups

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 sep-
arate branches on the right belong to the BI tree, which indicate two differences in topology with the ML tree.

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

Codon Count RSCU Codon Count RSCU
UUU(F) 199 1.66 UCU(S) 45 0.97
UUC(F) 41 0.34 UCC(S) 47 1.01
UUA\L) 267 2.66 UCA(S) 132 2.85
UUGIL) 29 0.29 UCG(S) 8 0.17
CUU(L) 110 1.09 CCU(P) 37 0.77
CUC(L) 32 G3? CCC(P) PF 0.57
CUA\(L) 143 1.42 CCA(P) 112 2530
CUG(L) 22 0.22 CCG(P) 15 0.31
AUU(!) 249 1 O8 ACU(T) 76 1.16
AUC() a2 0.45 ACC(T) 54 0.83
AUA(M) 194 1.67 ACA\T) 122 1.87
AUG(M) 38 0:33 ACG(T) 9 0.14
GUU(V) 53 1.16 GCU(A) 65 1.01
GUC(V) 22 0.48 GCC(A) 58 0.90
GUA(V) 88 1.93 GCA\A) Li? 84
GUG(V) 19 0.42 GCG(A) 18 0.28

*: stop codon.

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

Codon Count RSCU Codon Count RSCU
UAU(Y) 82 1.43 UGU(C) Ae 1.03
UAC(Y) 33 0.57 UGC(C) 14 0.97
UAA(*) 0 0.00 UGA(W) 91 1.70
UAG(*) 0 0.00 UGG(W) 16 0.30
CAU(H) 61 Laer CGU(R) ily 0.97
CAC(H) 35 0.73 CGC(R) 7 0.40
CAA(Q) 82 1.80 CGAI(R) 4l 2.34
CAG(Q) 9 0.20 CGG(R) 5 0.29
AAU(N) 85 1.12 AGU(S) 22 0.47
AAC(N) 67 0.88 AGC(S) 24 0.52
AAA\K) 78 LT AGA\(*) 0 0.00
AAG\(K) 10 0.23 AGG(*) 0 0.00
GAU(D) 43 130 GGU(G) ay 0.68
GAC(D) 23 0.70 GGC(G) 41 0.75
GAA\E) 74 1.54 GGA(G) 91 1.66
GAG(E) 22 0.46 GGG(G) 50 0.91

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

zse.pensoft.net

2/0 Zhang, Y. et al.: Mitogenome characteristics of P. flavomaculatus and biogeography of Hynobiidae

45.47
pe 22.27!
4.01
63.87 | 33.38

110.87

0.17 » Hynobius chinensis
2 Hynobius guabangshanensis
Hynobius maoershanensis
Hynobius yiwuensis
Hynobius amjiensis
Hynobius unisacculus
Hynobius quelpaertensis I
Hynobius yangi
25.65 Hynobius leechii

34.23 Hynobius nebulosus

Hynobius nigrescens
1.83 = Hynobius arisanensis
Hynobius formosanus
Hynobius kimurae
Batrachuperus pinchonii
Batrachuperus londongensis
Batrachuperus tibetanus
Batrachuperus yenyuanensis
Liua tsinpaensis
Liua shihi
Pseudohynobius flavomaculatus
Pseudohynobius jinfo
55.62 Pseudohynobius puxiongensis
Salamandrella tridactyla
Salamandrella keyserlingii
Pachyhynobius shangchengensis |
Batrachuperus mustersi |

= Batrachuperus gorganensis IV
45.27 Ranodon sibiri
162.79 a 0.39 » Ony, +
: Vv
181.17
ee 23.06 Echinotriton andersoni
Echinotriton chinhaiensis
10.06 Andrias davidianus Outgroups
Andrias japonicus
20.0
200 150 100 50 0 (Ma)

ee... Quaternary

Figure 4. Divergence time estimation of the Hynobiidae. Node bars indicate 95% Cls of the divergence time estimate. Chronostrati-
graphic scale data are derived from International Commission on Stratigraphy (https://stratigraphy.org/).

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 H 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 (1.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) hypothe-
sis, living hynobiids originated in northern China in the
Middle Cretaceous when the continent was relatively flat
and wet, and geological conditions facilitated dispersal

zse.pensoft.net

of hynobiids over broad geographical ranges. Since 50—
60 Ma, mountains such as the Himalayas and the Kun-
lun 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 geo-
graphical 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 emer-
gence of different species (Shen et al. 2020). Further-
more, speciation and geographical distribution of hynobi-
ids 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 w > 1, © = 0, and w = < 1 indicate that the effect of
natural selection on the change of PCGs is positive, neu-
tral, 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.

Zoosyst. Evol. 98 (2) 2022, 263-274 271
13PCGs
Ill
0.4500 II al
0.4000 SS IV
0.3500 I
0.3000 Vv
S 0.2500 ——
Z 0.2000
0.1500
0.1000
0.0500
0.0000 . ' ' : . ' . : 1 1 - : : - a 1 n : 1
a3 a5 RS o a5 ; Ro S eo my ne “s a ° 0 5 R A cs x & Ri ra ad
¢ oe oa ” a e e Rid es oe Po o s oe * & x se ye r re ro a” ee Oy of Ps eo se
. S P ey ye! SO vain! ae eat AD a) My Cae aT, a v eat x 0 \ 3
a \ eS < of Pago ¢ goo at we a9 8 » e ny Cae) er oFe <0 NN 5 &
S $ 3 6 NY RY) a. < OT oe eae 4 Rs) 3 a) a ay x >) we? rn ees &
wv Ss? * 0 ° \ S ¢ \ ~ of if RY « ny \ a yo So
y ay o : ) ay) \ <" : PN 5
x y e Ro AS es RY ws ws 3 ws oe Ro ° @ r # ot ye NY oY x »* xs a ‘ « F * i é Ra
Ri) wy s ey ~ oO “ ay af we oO wo RS < ae at we a Y RN x 3 x ¥ ny of x Ry
v ao ¥ er 6° gl gs . & et gl’ wt re to eh ge ot gr R & ey"
oO rg io Ro eg gl nS a ov OF
.) . * Ro s eo % a
ro g go ie)

Figure 5. The dN/dS ratios of the 13 concatenated PCGs for 32 Hynobiidae species. The Roman number represents the correspond-

ing partition in the phylogenetic tree.

Interestingly, one finding was surprising and concep-
tually 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 geo-
logical processes may have influenced rates of synony-
mous 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 1n 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. zhaoer-
mii and O. fischeri (0.39 Ma) diverged later, and their dN/
dS ratios are low. However, H. chinensis, H. guabang-
shanensis, H. maoershanensis, P. flavomaculatus, P. jinfo
do not appear to have speciated in a manner associated
with geologic events. Hynobius arisanensis & H. formo-
sanus may have arisen as a result of geological folds in
the central mountain range of Tatwan (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 lat-
er 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-va-
gility 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 substitu-
tions during the same period. Therefore, the dN/dS ra-
tios 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

zse.pensoft.net

272 Zhang, Y. et al.: Mitogenome characteristics of P. flavomaculatus and biogeography of Hynobiidae

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-di-
verging lineages of hynobiids were found to be higher
than those of later-diverging hynobiids. Together with
our biogeographic interpretation, we suggest that earli-
er 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 con-
ditions through time.

Conflicts of interest

The authors declare that they have no conflict of interest.

Acknowledgments

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 Commis-
sion (Grant No. KJQN202100503).

References

Altekar G, Dwarkadas S, Huelsenbeck JP, Ronquist F (2004) Parallel
metropolis coupled Markov chain Monte Carlo for Bayesian phy-
logenetic inference. Bioinformatics 20(3): 407-415. https://doi.
org/10.1093/bioinformatics/btg427

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.
https://doi.org/10.1038/290457a0

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. https://doi.
org/10.1093/genetics/108.1.237

Bernt M, Donath A, Juhling F, Externbrink F, Florentz C, Fritzsch G,
Putz J, Middendorf M, Stadler PF (2013) MITOS: Improved de
novo metazoan mitochondrial genome annotation. Molecular Phy-
logenetics and Evolution 69(2): 313-319. https://doi.org/10.1016/j.
ympev.2012.08.023

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.
https://do1.org/10.1038/s41559-018-0778-x

Boore JL (1999) Animal mitochondrial genomes. Nucleic Acids Re-
search 27(8): 1767-1780. https://doi.org/10.1093/nar/27.8.1767

Chai H, Du Y (2012) The complete mitochondrial genome of the pink
stem borer, Sesamia inferens, in comparison with four other Noc-
tuid moths. International Journal of Molecular Sciences 13(8):
10236-10256. https://doi.org/10.3390/ijms130810236

zse.pensoft.net

Crozier RH, Crozier YC (1993) The mitochondrial genome of the hon-
eybee Apis mellifera: Complete sequence and genome organization.
Genetics 133(1): 97-117. https://doi.org/10.1093/genetics/133.1.97

Dai L, Zhu B, Liu Q, Wei G, Liu C (2013) Characterization of the com-
plete mitochondrial genome of Bombyx mori strain H9 (Lepidop-
tera: Bombycidae). Gene 519(2): 326-334. https://doi.org/10.1016/j.
gene.2013.02.002

Darriba D, Taboada GL, Doallo R, Posada D (2012) jModelTest 2: More
models, new heuristics and parallel computing. Nature Methods
9(8): 772-772. https://doi.org/10.1038/nmeth.2109

Deng B (2016) Cenozoic mountain-building processes in the Daliang-
shan, 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) NOVOPIlasty: De novo
assembly of organelle genomes from whole genome data. Nucleic
Acids Research 45(4): e18. https://doi.org/10.1093/nar/gkw955

Edelaar P (2018) Ecological speciation: when and how variation among en-
vironments can drive population divergence. In: Bird Species, Spring-
er, Cham, 195—215. https://doi.org/10.1007/978-3-319-91689-7_ 11

Fang X (2017) Phased uplift of the Tibetan Plateau. Science & Technol-
ogy 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, Thal-
linger GG (2013) Complete mitochondrial DNA sequences of the
threadfin cichlid (Petrochromis trewavasae) and the blunthead
cichlid (7Jropheus moorii) and patterns of mitochondrial genome
evolution in cichlid fishes. PLoS ONE 8(6): 67048. https://doi.
org/10.1371/journal.pone.0067048

Gao K, Shubin NH (2003) Earliest known crown-group salamanders.
Nature 422(6930): 424-428. https://doi.org/10.1038/nature01491

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. https://doi.org/10.1002/ece3.5015

Huang H, Lin B, Guo C, Tang Z, Wu Z (2016) The complete mitochon-
drial genome of the Hynobius maoershanensis (Caudata: Hynobii-
dae). Mitochondrial DNA 27(1): 173-174. https://doi.org/10.3109/1
9401736.2013.878930

Kryazhimskiy S, Plotkin JB (2008) The population genetics of dN/dS.
PLOS Genetics 4(12): e1000304. https://doi.org/10.1371/journal.
pgen. 1000304

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.
https://doi.org/10.1101/gr.092759.109

Kumar S, Stecher G, Tamura K (2016) MEGA7: Molecular Evolution-
ary Genetics Analysis version 7.0 for bigger datasets. Molecular
Biology and Evolution 33(7): 1870-1874. https://doi.org/10.1093/
molbev/msw054

Kvist S, Utevsky S, Marrone F, Ben Ahmed R, Gajda L, Grosser C,
Huseynov M, Jueg U, Khomenko A, Oceguera-Figueroa A, PeSi¢
V, Pupins M, Rouag R, Salam N, Swiatek P, Trontelj P, Vecchioni
L, Muller C (2022) Extensive sampling sheds light on species-lev-
el diversity in Palearctic Placobdella (Annelida: Clitellata: Glossi-
phoniiformes). Hydrobiologia 849(5): 1-21. https://doi.org/10.1007/
$10750-021-04786-5

Zoosyst. Evol. 98 (2) 2022, 263-274

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.
https://doi.org/10.3109/19401736.2011.636437

Li Y, Wu M, Wang XL (2004) Phylogenetic relationships of Hynobii-
dae based on sequences of mitochondrial 16S ribosomal RNA gene.
Dong Wu Xue Bao 50(3): 464-469.

Li J, Fu C, Lei G (2011) Biogeographical Consequences of Cenozo-
ic Tectonic Events within East Asian Margins: A Case Study of
Hynobius Biogeography. PLoS ONE 6(6): e21506. https://doi.
org/10.1371/journal.pone.0021506

Liu JQ, Chen SS, Guo WF, Sun CQ, Zhang ML, Guo ZF (2015)
Research advances in the Mt. Changbai volcano. Bulletin of Miner-
alogy. Petrology and Geochemistry 34(4): 710-723.

Mccormack JE, Huang H, Knowles LL (2009) Maximum Likelihood Es-
timates of Species Trees: How Accuracy of Phylogenetic Inference
Depends upon the Divergence History and Sampling Design. System-
atic Biology 58(5): 501-508. https://doi.org/10.1093/sysbio/syp045

Meng QR (2017) Origin of the Qin ling Mountains. Scientia Sinica Ter-
rae 47: 412-420. https://doi.org/10.1360/N072016-00422

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. https://doi.org/10.1007/BF02 100203

Peng R, Zhang P, Xiong J, Gu H, Zeng X, Zou F (2010) Rediscov-
ery of Protohynobius puxiongensis (Caudata: Hynobiidae) and its
phylogenetic position based on complete mitochondrial genomes.
Molecular Phylogenetics and Evolution 56(1): 252-258. https://doi.
org/10.1016/j.ympev.2009.12.011

Perna NT, Kocher T (1995) Patterns of nucleotide composition at fourfold
degenerate sites of animal mitochondrial genomes. Journal of Molec-
ular Evolution 41(3): 353-358. https://doi.org/10.1007/BF01215182

Rambaut A, Drummond AJ, Xie D, Baele G, Marc A (2018) Posterior
Summarization in Bayesian Phylogenetics Using Tracer 1.7. System-
atic Biology 67(5): 901-904. https://do1.org/10.1093/sysbio/syy032

Reilly SM (1987) Paradactylodon: A Junior Synonym for Ba-
trachuperus. Amphibia-Reptilia 8(3): 283-284. _ https://doi.
org/10.1163/156853887X00324

Ronquist F, Huelsenbeck JP (2003) MrBayes 3: Bayesian phylogenetic
inference under mixed models. Bioinformatics 19(12): 1572-1574.
https://doi.org/10.1093/bioinformatics/btg 180

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. https://doi.
org/10.1007/s11427-012-4348-1

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 Evo-
lution 7(6): 1869-1881. https://doi.org/10.1002/ece3.2737

Shen Y, Yang N, Liu Z, Chen Q, Li Y (2020) Phylogenetic perspec-
tive on the relationships and evolutionary history of the Acipens-
eriformes. Genomics 112(5): 3511-3517. https://do1.org/10.1016/).
ygeno.2020.02.017

Stock M, Fakharzadeh F, Kuhl H, Rozenblutkoscisty 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,

213

and Transcriptomic Data. Genes 10(4): 306. https://doi.org/10.3390/
genes10040306

Su T, Liang A (2019) Comparative analysis of seven mitochondrial ge-
nomes of Phymatostetha (Hemiptera: Cercopidae) and phylogenetic
implications. International Journal of Biological Macromolecules
125: 1112-1117. https://doi.org/10.1016/).ijbiomac.2018.12.174

Tang D, Xu T, Sun Y (2015) Hynobiidae origin in middle Cretaceous
corroborated by the new mitochondrial genome of Hynobius chin-
ensis. Marine Genomics 22: 37—44. https://doi.org/10.1016/j.mar-
gen.2015.03.006

Thorn R, Raffaélli J (2001) Les Salamandres de l’Ancien Monde. Paris:
Boubée, 380 pp.

Uddin A, Chakraborty S (2019) Codon Usage Pattern of Genes In-
volved in Central Nervous System. Molecular Neurobiology 56(3):
1737-1748. https://doi.org/10.1007/s12035-018-1173-y

Wang C, Dai J, Zhao X, Li Y, Graham SA, He D, Ran B, Meng J (2014)
Outward-growth of the Tibetan Plateau during the Cenozoic: A review.
Tectonophysics 621: 1-43. https://doi.org/10.1016/j.tecto.2014.01.036

Wolstenholme DR (1992) Animal mitochondrial DNA: structure and
evolution. International Review of Cytology-a Survey of Cell Biol-
ogy 141: 173-216. https://doi.org/10.1016/S0074-7696(08)62066-5

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 Long-
maxi formation in east Sichuan Wuling mountain area. China Uni-
versity of Geosciences.

Xiong JL, Sun P, Zhu WW, Liu XY (2009) Review on the Taxonomic
Variation of Genus and Species in Chinese Hynobiids. Sichuan Jour-
nal 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 Technolo-
gy 11(1): 55-58.

Xu J (2002) Some taxonomic problems of the Hynoiidae. Acta Scien-
tiarum 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. Sheng-
taixue 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 systemat-
ics of the hynobiid salamanders of the Liua—Pseudohynobius complex:
Molecular and phylogenetic perspectives. Biochemical Systematics and
Ecology 34(6): 467-477. https://doi.org/10.1016/j.bse.2006.01.006

Zhang Y (2015) Study of tectonic thermal evolution history and up-
lift 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. https://doi.
org/10.1016/S1055-7903(03)00059-9

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.
https://doi.org/10.1073/pnas.0602325103

zse.pensoft.net

274 Zhang, Y. et al.: Mitogenome characteristics of P. flavomaculatus and biogeography of Hynobiidae

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-Neo-
gene stratigraphic realm and sedimentary sequence of the Qing-
hai-Tibet Plateau and their response to uplift of the plateau. Science
China-earth Sciences 53(9): 1271-1294. https://doi.org/10.1007/
$11430-010-4048-2

Zhang Z, Dai W, Dai D (2013) Synonymous codon usage in TTSuV2:
Analysis and comparison with TTSuV1. PLoS ONE 8(11): e81469.
https://doi.org/10.1371/journal.pone.008 1469

Zhao EM, Hu QX (1984) Studies on Chinese Tailed Amphibians. Sich-
uan 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. https://doi.
org/10.1371/journal.pone.0078064

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. https://doi.org/10.1093/molbev/msr072

Zhou Z, Xie F, Jiang JP, Zheng ZH (2004) Research advances of amphib-
ian declines. Chinese Journal of Applied and Environmental Biolo-
gy 10(1): 128-132. https://doi.org/10.3724/SP.J.1145.2012.00128

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. https://do1.org/10.1002/ece3.8764

Supplementary material |

Figure SI

Authors: 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 la-
belled with the abbreviations of their corresponding ami-
no acids. The lines between the English letters indicate
the phosphodiester bonds and Watson—Crick pairing.

Copyright notice: This dataset is made available under
the Open Database License (http://opendatacommons.
org/licenses/odbl/1.0/). The Open Database License
(ODbL) is a license agreement intended to allow us-
ers 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.

Link: https://do1.org/10.3897/zse.98.66578.suppl 1

zse.pensoft.net

Supplementary material 2
Figure S2

Authors: 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.

Copyright notice: This dataset is made available under
the Open Database License (http://opendatacommons.
org/licenses/odbl/1.0/). The Open Database License
(ODbL) is a license agreement intended to allow us-
ers 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.

Link: https://doi.org/10.3897/zse.98.66578.suppl2

Supplementary material 3

Tables S1, S2

Authors: 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. Composi-
tion and skewness of the mitogenome for 36 Caudata
species.

Copyright notice: This dataset is made available under
the Open Database License (http://opendatacommons.
org/licenses/odbl/1.0/). The Open Database License
(ODbL) is a license agreement intended to allow us-
ers 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.

Link: https://doi.org/10.3897/zse.98.66578.suppl3