Complete Mitogenome of Oreolalax omeimontis Reveals Phylogenetic Status and Novel Gene Arrangement of Archaeobatrachia

Species of the genus Oreolalax displayed crucial morphological characteristics of vertebrates transitioning from aquatic to terrestrial habitats; thus, they can be regarded as a representative vertebrate genus for this landing phenomenon. But the present phylogenetic status of Oreolalax omeimontis has been controversial with morphological and molecular approaches, and specific gene rearrangements were discovered in all six published Oreolalax mitogenomes, which are rarely observed in Archaeobatrachia. Therefore, this study determined the complete mitogenome of O. omeimontis with the aim of identifying its precise phylogenetic position and novel gene arrangement in Archaeobatrachia. Phylogenetic analysis with Bayesian inference and maximum likelihood indicates O. omeimontis is a sister group to O. lichuanensis, which is consistent with previous phylogenetic analysis based on morphological characteristics, but contrasts with other studies using multiple gene fragments. Moreover, although the duplication of trnM occurred in all seven Oreolalax species, the translocation of trnQ and trnM occurred differently in O. omeimontis to the other six, and this unique rearrangement would happen after the speciation of O. omeimontis. In general, this study sheds new light on the phylogenetic relationships and gene rearrangements of Archaeobatrachia.


Introduction
The mitochondrial genome (mitogenome) is an important model system for studying molecular evolution, phylogeny, and genome structure [1,2]. Most vertebrate mitogenomes are characterized by several features shared with most eukaryote mitogenomes, such as a circular, double-stranded DNA molecules, which generally range from 15 to 20 kb [3]. This smaller organelle genome consists of 13 protein-coding genes (PCGs), two ribosomal RNA (rRNA) genes, 22 transfer RNA (tRNA) genes, and an AT-rich region, known as the control region (CR), which is involved in the initiation of transcription and replication of the mitogenome [4]. Due to their small size, cellular abundance, typical maternal inheritance, compact gene arrangement, and high rate of evolution, mitogenomes have been extensively used to test hypotheses about microevolution, to study population structure, phylogeny, and phylogeography at various taxonomic levels, and to identify cryptic species, which complements new developments using nuclear genes [1,5,6].
The genus Oreolalax (Anura, Megophryidae), which was reevaluated by Fei et al. [7] and has been accepted [8], currently consists of 19 species and is believed to originate from southwestern China [8] but has been found unexpectedly in Vietnam [9]. Wei et al. [10] considered that the genus Oreolalax contains the crucial morphological characteristics of species transitioning from aquatic to terrestrial habitats. Therefore, Oreolalax has become an important model for answering basic questions of morphology, development, and biogeographic evolution [11].
Six Oreolalax mitogenomes have been published [12][13][14][15][16] and have been shallowly analyzed for their phylogenetic relationships, which has led to a preliminary understanding of the phylogenetic relationships of Oreolalax. Furthermore, although gene rearrangement is rare in non-neobatrachians (i.e., archaeobatrachians) [17][18][19], specific gene rearrangements were discovered in all six Oreolalax mitogenomes. As for O. omeimontis (Omei toothed toad [20] or Omei lazy toad [21], common name), which can be distinguished from other Oreolalax species by its characteristic triangular markings between the eyes, this native Chinese species, mainly distributed in the Emei Mountains [22], somehow has only a few genetic fragments but no complete mitogenomic data. This further hinders the study of the evolution of gene rearrangements in the mitogenome of the genus Oreolalax and the Archaeobatrachia suborder.
The present phylogenetic status of O. omeimontis is still controversial and needs to be addressed, as morphology-based phylogenetic studies suggest that O. omeimontis and O. lichuanensis are sister branches [10], while phylogenetic studies based on molecular data (gene fragments) suggest that O. omeimontis was a sister to O. multipunctatus [9,23,24]. However, morphological characteristics change with the development of an individual, which poses methodological limitations [25], and differences in morphological characteristics may also be accompanied by variation in mitogenomes [26,27]. Meanwhile, previous research has also suggested that longer DNA sequences, such as a complete mitogenome, can provide more phylogenetic information than shorter gene fragments, which may more accurately reveal phylogenetic relationships and systematic taxonomic status [28][29][30]. Therefore, further research into the mitogenomic features, gene rearrangements and phylogenetics of O. omeimontis was carried out here to explore this poorly understood species and supplements the existing molecular data for Oreolalax taxa.
To fill the gap in our understanding of O. omeimontis and provide deep-level phylogenetic analysis for Oreolalax, this study sequenced the complete mitogenome of O. omeimontis, and described its characteristics and gene rearrangement. Moreover, we focused on the comparison between the seven mitogenomes of Oreolalax species to explore their features; maximum likelihood (ML) and Bayesian inference (BI) methods were employed to perform phylogenetic analyses to provide further information on the genus Oreolalax. This provides a valuable resource for further studies on the genetic diversity and phylogenetic analysis of Oreolalax.

Sample Collection and DNA Extraction
O. omeimontis specimens were collected from Wawushan National Forest Park, Sichuan, China in 2015. Tissue samples were preserved in 99% ethanol until DNA extraction. Total genomic DNA was extracted from muscle tissue according to protocols of the Ezup Column Tissue Genomic DNA Purification Kit (Sangon Biotech, Shanghai, China).

PCR Amplification and Sequencing
Twelve paired primers were used to amplify the complete mitogenome of O. omeimontis according to Zhang et al. [31]. PCR amplifications were performed with 2 µL of genomic DNA, 1 µL of each primer (10 µmol/L), 6 µL of ddH 2 O, 10 µL of 2×Taq Master-Mix Kit (Sangon Biotech, Shanghai, China). PCR reactions were performed under the following conditions: 5 min at 94 • C, followed by 35 cycles of 40 s at 94 • C, 40 s at 48 to 55 • C (depending on the primer combination), elongation at 72 • C for 1 to 3 min (depending on the fragment length), and a final extension at 72 • C for 10 min. PCR products were electrophoresed on a 1.0% agarose gel, then sequenced by Tsingke Biotechnology (Chengdu, China). DNA sequencing was performed on an ABI 3730 automatic DNA sequencer (Applied Biosystems, Foster City, CA, USA).

Sequence Assembly, Annotation, and Analysis
Each sequence was checked using the NCBI BLAST database (https://blast.ncbi. nlm.nih.gov/Blast.cgi, accessed on 10 May 2020), and then was assembled using both MEGA 7 [32] and SeqMan programs (Lasergene v7.1.0; DNAStar, Inc., Madison, WI, USA). A total of thirteen protein coding genes (PCGs) were first identified using the NCBI ORF Finder (https://www.ncbi.nlm.nih.gov/orffinder/, accessed on 10 May 2020) to specify the vertebrate mitochondrial genes. PCG, rRNA, and tRNA genes were annotated by the MITOS web server (http://mitos2.bioinf.uni-leipzig.de/index.py, accessed on 10 May 2020) [33] and tRNAscan-SE (http://trna.ucsc.edu/tRNAscan-SE/, accessed on 10 May 2020) [34] was used to compare and ensure the accuracy of the results; furthermore, both MITOS and tRNAscan-SE was used for recognition and prediction of tRNA gene secondary structure. The circular graphical map of the mitogenome was drawn using CGview web server (https://cgview.ca/, accessed on 10 May 2020) [35]. The relative synonymous codon usage (RSCU), start/stop codon, codon usages, and composition of nucleotides were analyzed in MEGA 7. AT-skew and GC-skew were calculated using the following formulas: AT-skew = (A -T) / (A + T) and GC-skew = (G -C) / (G + C) [36]. In addition, the same methods were used to reanalyze the published mitogenomes of six other Oreolalax species to supplement previously unanalyzed data.

Phylogenetic Analysis
To determine the phylogenetic position of O. omeimontis, 31 frogs (including O. omeimontis of this study) representing eight Anura families and three outgroup species (Table S1) were included for this analysis. The complete sequences of 13 PCGs of all species with deletions in the stop codons were aligned through MAFFT v7.313 [37] using default parameters. These 13 single-gene sequence alignments were concatenated using MEGA 7 and then used for downstream phylogenetic purposes. PartitionFinder v2.1.1 [38] was used to select the nucleotide substitution model of 13PCGs for the construction of the BI phylogenetic tree, and ModelFinder [39] was used to select for the ML phylogenetic tree, respectively. For ML and BI analysis, the nucleotides of 13 PCGs were codon partitioned and the best-fit model of each partition was selected (Table S2). Phylogenetic analysis was performed with the ML and Bayesian inference (BI) methods using IQ-TREE [40] and Mrbayes v.3.2.2 [41], respectively. For ML analysis, 1000 bootstrap replicates were used to calculate the bootstrap of the program. In the BI analysis, two independent runs for 10 7 Markov chain Monte Carlo (MCMC) generations with four chains and sampling trees occurred every 1000 generations. The first 25% of trees were discarded as burn-in samples and the remaining trees were used to generate Bayesian consensus trees. The above software was running on the integrated and scalable desktop platform of PhyloSuite v1.1.16 [42]. The iTOL v4 website (https://itol.embl.de/, accessed on 10 May 2020) [43] was used to visualize and edit the results of both ML and BI trees.

Mitogenome Organization
The complete mitogenome of O. omeimontis was 17,266 bp in length (GenBank accession numbers: OP722573). The gene content was typical of other Oreolalax mitogenomes, which contained 13 PCGs, 2 rRNA genes, 23 tRNA genes (including an additional copy of trnM), a non-coding region known as the CR, and three intergenic regions ( Figure 1). Eight tRNA genes and nad6 were located on the L-strand (-), and the remaining genes (including the additional copy of trnM) were encoded by the H-strand (+) ( Table 1). The overall base composition of O. omeimontis was A: 27.6%, T: 31.3%, C: 26.2%, and G: 14.8%. The nucleotide composition of each part of the O. omeimontis mitogenome was A + T rich, ranging from 58.2% to 63.9% for the whole mitogenome of 13 PCGs, 2 rRNAs, 23 tRNAs and CR. The AT skewness (−0.063) and GC skewness (−0.278) values were both negative, indicating bias towards Ts rather than As, and bias towards Cs rather than Gs for the whole mitogenome of O. omeimontis. Table 2 shows the base composition of seven species of Oreolalax according to their mitogenomes, PCGs, tRNA genes, rRNA genes and control region. Compared to the six other Oreolalax species, O. omeimontis had the lowest A + T content in the complete mitogenome. Moreover, O. omeimontis had the lowest GC-skew values in tRNAs and PCGs, indicating that O. omeimontis had the lowest Gs content in tRNAs and PCGs out of the seven species of Oreolalax. In addition, the seven Oreolalax species' AT-skew and GC-skew data had the same positive and negative values in tRNAs, PCGs, rRNAs, and CR, indicating that they were roughly the same in base biasing.

Protein-Coding Genes and Codon Usage
The total length of 13 PCGs for O. omeimontis was 11,389 bp. The coding regions ranged in size from 168 (atp8) to 1821 bp (nad5), accounting for approximately 66% of the entire mitogenome. Despite nad6 being located on the L-strand (-), the remainder of the genes (nad1, nad2, nad3, nad4l, nad4, cox1, cox2, cox3, atp8, atp6, nad5, and cob) are all located on the H-strand (+). For the seven Oreolalax species, most 11 PCGs started with an ATN codon (ATG for nad1, nad2, cox2, atp8, atp6, cox3, nad4l, nad4, nad5, and cob; ATA for nad3) (Table S3). In contrast, most cox1 and nad6 started with the GTG codon. Ten PGCs ended with the stop codon TAA or TAG and included an incomplete stop codon T-or TA-, while cox1, nad5 (except for O. jingdongensis) and nad6 terminated with the complete stop codon AGG or AGA (Table S3). These were all canonical mitochondrial initiation and termination codons for vertebrate mitogenomes [44,45]. The relative synonymous codon usage (RSCU) values of the O. omeimontis mitogenome are shown in Figure 2 and Table (Table S4). Obviously, overrepresentation of the Us (Ts) in frequently used codons indicates that the codon usage of PCGs contributed to both the A + T bias, and that Ts (Us) are the most abundant base of Oreolalax mitogenomes (

Transfer RNA
A total of 23 tRNA genes were found in the mitogenome of O. omeimontis, including an additional copy of trnM, which is consistent with other studies of Oreolalax species [13][14][15]. The tRNA genes ranged in size from 64 (trnC) to 75 bp (trnL2). Among them, eight tRNA genes were located on the L-strand (-) and the remaining 15 are located on the H-strand (+). Most of the tRNA genes were found to fold into the typical clover-leaf structure (Figure 3). We found a high degree of similarity between two trnM of O. omeimontis and between all the trnM form the other six Oreolalax species. In addition, the secondary structures of the two trnM sequences were predicted successfully. They had the same anticodon of CAU and could fold into the typical clover-leaf structure (Figure 3). This result indicates that neither of the two trnM sequences had degenerated, although a small number of non-complementary and U-G base pairs existed in the stem regions. Stem mismatches are common within mitochondrial tRNA genes and are repaired via post-transcriptional editing [46,47]. Furthermore, a lack of the dihydrouridine (DHU) in the arm of trnS1 does not seem to affect its function, according to previous studies [48,49].

Ribosomal RNA and Control Region
For the O. omeimontis, rrnS was 935 bp long, which was located between trnF and trnV; rrnL was 1600 bp long and was present between trnV and trnL2. For rRNAs, the total A + T content was 58.2%, with a positive AT-skew (0.124) and a negative GC-skew (−0.105) ( Table 2), indicating that the content of As is highest and the content of Gs is lowest. In contrast, O. omeimontis had the lowest A + T content (58.2%) among the Oreolalax species (Table 2). Consistent with the most metazoan mitogenomes (Table 1), O. omeimontis contains only one control region, which was located between the trnW and trnF genes with 1166 bp in size. The CR had the highest A + T content (63.9%) in the mitogenome, with both negative AT-skew (−0.095) and GC-skew (−0.291) ( Table 2). Moreover, in all seven Oreolalax species, the CR had the highest A + T content in mitogenomes with both negative AT-skew and GC-skew, indicating that bias exists toward Ts and Cs in Oreolalax in CRs.

Phylogenetic Analysis
Phylogenetic analyses were carried out based on the concatenated alignment of nucleotides sequences of 13 PCGs covering eight families and 31 species of Archaeobatrachia; three species were used as outgroups (Table S1). The final concatenated alignment of our mitogenome dataset for 34 species contained 11,478 nucleotide positions. All the topological structures from the Bayesian and ML analyses were consistent (Figure 4). The present phylogenetic analysis of archaeobatrachians shows ((((((((Megophryidae + Pelobatidae) + Pelodytidae) + Scaphiopodidae)) + Pipidae)) + (Bombinatoridae + Alytidae)) + Leiopelmatidae) (Figure 4), which is consistent with previous studies [31,50]. The six genera of Megophryidae were clustered as (((Oreolalax + Leptobrachium) + Scutiger) + ((Xenophrys + Megophrys) + Atympanophrys)), which indicates that Oreolalax is sister group to Leptobrachium, and forms a branch with Scutiger ( Figure 4). This result is consistent with Feng et al. [50] but inconsistent with Pyron et al. [25], who proposed that Oreolalax should be a sister group to Scutiger but not Leptobrachium. We believe that the incongruence between these two studies is that Feng et al. [50] used a longer sequence alignment with more genes (around 88,000 nt of aligned sequences from 95 nuclear protein-coding genes), while Pyron et al. [25], which predates Feng et al. [50] by six years, used only 12 genes (nine nuclear plus three mitochondrial genes in 12,712 nt) for phylogenetic purposes. As a result, both the previous study and ours were able to assemble more divergent loci to generate a more accurate phylogenetic structure.
Seven  (Figure 4). O. rhodostigmatus was placed at the basal position of the genus Oreolalax, which is consistent with previous studies that employed multi-markers [9,24,25]. O. major was a sister to O. xiangchengensis and O. omeimontis forms a sister relationship with O. lichuanensis, which was well supported by statistical data (bootstrap value, BSV = 100; Bayesian posterior probability, BPP = 1). Previous studies based on 13 mitochondrial PCGs [13][14][15] or partial nuclear and mitochondrial [24,25] data had the consistent view that O. major and O. xiangchengensis form a sister group. However, previous phylogenetic analyses based on multiple nuclear and mitochondrial gene fragments reconstructed O. omeimontis and O. multipunctatus as a sister branch [9,24,25], whereas our results based on 13 PCGs from the complete mitogenome were consistent with the phylogenetic results based on 21 morphological characters [23], i.e., O. omeimontis and O. lichuanensis were reconstructed as sister branch. Similar to the previous studies suggesting the use of longer sequences to obtain more accurate phylogenetic results [29], we believe that our study using 13 mitochondrial PCGs provided a more reliable phylogenetic relationship than research using a small number of gene fragments.
In general, our results demonstrated a different view to the studies with only a small number of gene fragments, as well as supported traditional approaches of using morphological data for phylogenic purposes, and also complemented phylogenetic analyses involving genus Oreolalax. Thus, more taxon sampling and more sequence data, including both nuclear and mitochondrial, as well as morphological markers, are required to obtain a more robust phylogenetic framework for Oreolalax taxa.
In this study, the gene rearrangement events among 31 archaeobatrachians mitogenomes revealed five gene arrangement types (type I to type V); only the O. omeimontis mitogenome was recognized as type V ( Figure 5). All species except Leiopelmatidae (type I) and partial Megophryidae species (types III, IV and V) have the typical (traditional) vertebrate mitogenome gene order (type II). Compared to typical mitogenome organization (type II), in the basal lineage of Leiopelmatidae, the two adjacent genes of nad6 and trnE has been translocated from the canonical location (type II) between the nad5 and cob genes to a region downstream of trnT (type I). Moreover, the translocation of trnI, trnW, and trnF seems to have occurred in Scutiger ningshanensis (type III), which is unique for other archaeobatrachians ( Figure 5). Additionally, duplication of trnM and translocation of trnW had both occurred in six Oreolalax species and two Leptobrachium species (type IV), and both Oreolalax and Leptobrachium were closely related in phylogenetic relationships, these gene rearrangement events may have occurred to the common ancestor of these two genera. Moreover, the translocations of trnQ and trnM in O. omeimontis seem to result from an independent event, as O. omeimontis's arrangement (type V) is unique to other six Oreolalax species ( Figure 5) Gene rearrangements have been shown to occur independently in multiple lineages, which are relatively rare and random [19,57,58]. To explore whether there are more rearrangement events in the genus Oreolalax as well as Archaeobatrachia suborder, and to further investigate the mechanism of mitogenome rearrangement and evolution of Anura, further studies and more species of Oreolalax should be sequenced. Such follow-up studies could lead to new discoveries and insights into the evolution of mitogenome reorganization.

Conclusions
This study has presented the characterization of the complete mitogenome of O. omeimontis, which complements the molecular data of Oreolalax taxa. In addition, the gene arrangement of O. omeimontis, surprisingly, was revealed to be different from the other six published Oreolalax species. Furthermore, phylogenetic analyses were conducted for some Archaeobatrachia species and confirmed the phylogenic position of O. omeimontis as a sister species to O. lichuanensis by using 13 mitochondrial PCGs. In general, this study provides a novel insight into the mitochondrial gene arrangements and phylogenetic framework of the Archaeobatrachia suborder, and also facilitates future studies on the evolution of vertebrates. Meanwhile, further sampling from different taxonomic grades for mitogenomic and nuclear data will be helpful in understanding the phylogenetic and evolutionary implications among the species of Archaeobatrachia.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/genes13112089/s1, Table S1: List of the mitogenomes of species analyzed in this study and their GenBank accession numbers; Table S2: Best partitioning schemes selected by ModelFinder (13PCGS -ML) and PartitionFinder2 (13PCGS -BI), respectively; Table S3. Usage of start and stop codons in the mitogenomes of seven species of Oreolalax; Table S4. Condon number and RSCU of seven Oreolalax species mitochondrial PCGs. The asterisk in parentheses indicates the stop codon.