Three New Ranidae Mitogenomes and the Evolution of Mitochondrial Gene Rearrangements among Ranidae Species

Jiandong YANG, Jiaojiao YU, Jiabin LIU, Ming ZHOU, Biao LI and Bo OUYANG

1 College of Animal Science and Technology, Sichuan Agricultural University, Chengdu 611130, Sichuan, China

2 Sichuan Key Laboratory of Conservation Biology for Endangered Wildlife, Chengdu Research Base of Giant Panda Breeding, Chengdu 610081, Sichuan, China

1. Introduction

Previous studies had revealed that the gene organization in vertebrate mitogenomes is conserved and that the mitochondrial D-loop region and the 37 genes were arranged in same manner among vertebrates (Andersonet al., 1981; Roeet al., 1985; Tzenget al., 1992; Zardoyaet al., 1995). However, numerous gene rearrangements in the mitogenome can independently evolve (Alamet al.,2010; Chenet al., 2011; Desjardins and Morais, 1990;Kurabayashiet al., 2006, 2008, 2010; Liuet al., 2005;Mindellet al., 1998; Moritz and Brown, 1987; Sanoet al., 2005; Suet al., 2007; Zhanget al., 2013). Gene rearrangements involve duplications, losses, translocation,inversion, and/or shuffling of the D-loop region (also known as the control region), the replication origin of the light strand (OL) and the codon genes (including rRNA genes, tRNA genes and protein-coding genes).Although distinct mitogenome structural features have been reported for some amphibians, most amphibians(including caecilians, salamanders, archaeobatrachians,and mesobatrachians) generally conform to the typical Vertebrate-type mitochondrial gene arrangement (Liuet al., 2016; Mueller and Boore, 2005; Pabijanet al.,2008; San Mauroet al., 2004, 2006, 2014; Xiaet al.,2010; Zhanget al., 2008; Zhang and Wake, 2009).Surprisingly, the gene arrangements in the neobatrachian group are especially diverse and complex, and notably,their four tRNA genes (LTPF-trn) are commonly rearranged, which is distinguishable from the vertebrate ancestral gene order (Kurabayashiet al., 2010; Sumidaet al., 2001; Xiaet al., 2014).

The vertebrate mitochondrial rearrangements appear to be unique, random, generally rare events that are exceptionally unlikely to arise independently in independent evolutionary lineages (Boore and Brown, 1998; Liu and Huang, 2010; Xiaet al., 2010,2014), although a few convergent or parallel gene rearrangements have been observed in vertebrate mtDNAs(Morrisonet al., 2002; Weiet al., 2014). The exceptional mitochondrial gene rearrangement has been thought to have significant implication for animal phylogenetic analysis and is considered a powerful phylogenetic marker also applicable to explore phylogenetic relationships among various groups at different taxonomic levels (Boore and Brown, 1998; Maceyet al., 1997; San Mauroet al., 2004, 2014; Weiet al., 2014; Xiaet al.,2010; Xueet al., 2016; Zhanget al., 2008, 2009, 2013).For example,Odorrana tormota, a species famous for its ultrasonic communication, was previously regarded as a member ofAmolops(Frost, 2017). However, this frog shares the same mitochondrial gene arrangement (thetrnHwas translocated to D-loop downstream, forming aHLTPF-trncluster) with mostOdorranafrogs, not with theAmolopsfrogs (conventionalLTPF-trncluster) (Suet al., 2007).

The family Ranidae, also known as ranid frogs, is one of the most species-rich and fascinating groups of vertebrates (Cheet al., 2007; Liet al., 2014; Frost,2017). Ranidae represents one of the main components of Neobatrachia and contains approximately 380 described species, belonging to 23–24 genera (AmphibiaWeb, 2017;Frost, 2017). A total of 31 complete and 12 near-complete ranid mitochondrial genomes have been submitted to GenBank, and many novel gene rearrangement types have also been discovered (e.g. Liet al., 2014, 2016a, b;Kurabayashiet al., 2010; Suet al., 2007). Kurabayashiet al.(2010) reported the partial or complete mtDNAs of 10 ranids and found most mitogenomes were different from the typical Neobatrachian-type gene arrangement.The diversity of mitochondrial gene arrangements in ranid species is unexpected high (Kurabayashiet al., 2010).

Here, we decode the mitochondrial genomes of three ranid frogs, conduct comparative genome analysis with all available Ranidae mitogenome sequences submitted to GenBank, and perform the phylogenetic analysis among Ranidae species. Our aim was to conduct an in-depth investigation, including examining the phylogenetic relationships, redescribing the novel mitogenome structures, analyzing exhaustively the genome reorganization types, and inferring the possible mechanisms and evolutionary pathways of gene rearrangements as well as its systematic implication among ranid frogs. Our study helps to understand mitogenome evolution and phylogenetic relationships of Ranidae species.

2. Materials and Methods

2.1. Specimen collection, DNA extraction, and PCR amplificationSpecimens ofRana kukunoris,R. chaochiaoensisandR. omeimontiswere obtained from Zoige County (33.57066° N, 102.96348° E, 3446 m a.s.l.), Shimian County (29.02461° N, 102.38626°E, 2 085 m a.s.l.), and Yucheng District (29.97900° N,102.98117° E, 618 m a.s.l.) in Sichuan Province, China,respectively, and stored at -80°C. A TaKaRa MiniBEST Universal Genomic DNA Extraction Kit Ver.5.0 (Takara,Dalian, China) was used to extract total genomic DNA from a frozen tissue sample of the thigh muscle according to the detailed manufacturer’s protocol. Primer sets used to amplify the entire mitogenomes of the threeRanaspecies are shown in Table S1.

2.2. Sequence assembly and annotationThe overlapping sequence fragments were assembled by the program Seqmen (DNAstar, Madison, WI, USA).The annotations of rRNA genes (rRNAs), tRNAs,protein coding genes (PCGs) and D-loop region and the definitions of their respective gene boundaries were performed by the MitoAnnotator service (http://mitofish.aori.u-tokyo.ac.jp/annotation/input.html). The ARWEN program (http://mbio-serv2.mbioekol.lu.se/ARWEN/) was also utilized to infer the tRNAs via their proposed cloverleaf secondary structure and anticodon sequences. All annotation results were verified via alignment with homologous regions from other reportedRanamitochondrial genomes. Finally, the mitochondrial genetic diagrams were generated by the OGDRAW program (http://ogdraw.mpimp-golm.mpg.de).

2.3. Data collectionWe downloaded 32 complete and 12 partial Ranidae mitochondrial genomes from GenBank(Table 1). Eight non-Ranidae mitogenomes were used as out-groups in the phylogenetic analysis. The taxonomic names of all species were based on ‘Amphibian Species of the World 6.0’ (Frost, 2017). There were many errors in some mitogenome annotations previously submitted to GenBank, and these mitogenome sequences should be re-annotated in systematic or comparative research(Cameron, 2014). In order to avoid interference caused by these errors in our subsequent analysis, we reanalyzed all sequences using the online services MitoAnnotator and ARWEN. The important corrections were listed in Figure S1.

2.4. Genome rearrangement analysisWe compared and analyzed re-annotated mitogenomes, together with the three newRanafrog data, with respect to mitogenome gene order (Chenet al., 2011). The definition of mitogenome organization types is based on the comparative results. To clarify, if the gene arrangements of the new mitogenome deviate from the typical Vertebrate-type gene arrangement (Type A) and the typical Neobatrachian-type gene arrangement (Type B), we will divide it into a new type (Figure 1). The long intergenic spacer frequently found in the closely related species and the pseudogene are also taken into account.If we cannot determine that the long intergenic spacer(more than 20 bp in size) frequently found in the closely related species is a pseudogene via homologous sequence alignments, for convenience, we will temporarily call it as“gap” in this study.

2.5. Phylogenetic tree analysisFirstly, all termination codons of 13 PCGs nucleotide sequences were manually deleted. Then, the remaining fragments of each PCG were separately aligned based on their translated amino acid sequences by Muscle implemented in MEGA6.06(Tamuraet al., 2013), and the two rRNAs sequences were separately aligned by ClustalX2 (Larkinet al.,2007). Subsequently, all ambiguous alignment regions were trimmed by the Gblocks Server (http://molevol.cmima.csic.es/castresana/Gblocks_server.html), the type of sequence was set to Codons (for PCGs) or DNA (for rRNAs) and all options for a less stringent selection were selected. Finally, the 15 trimmed alignments were concatenated into a single dataset to infer the phylogenetic relationships of Ranidae. For the concatenated sequence matrix, two phylogenetic trees were constructed using both Bayesian inference (BI) and maximum likelihood(ML) approaches. The ML analysis was conducted by PhyML3.1 (Guindonet al., 2010) under the GTR + I +G evolutionary model determined by jModelTest2.1.5(Darribaet al., 2012), with 100 replicates for the nonparametric bootstrap analysis. The BI analysis was performed by MrBayes3.2.2 (Ronquist and Huelsenbeck,2003). For the BI analysis, we firstly partitioned the data into 15 partitions by gene, and then used jModelTest2.1.5 to select the best- fit model of nucleotide substitution for each partition with the Bayesian Information Criterion,which was preferred for model selection (Luoet al.,2010). We performed two independent runs for 5 000 000 generations, sampled every 1 000 generations,conservatively discarded the first 25% of generations as burn-in, and visualized the majority-rule (＞50%)consensus trees using FigTree1.4.2 (http://tree.bio.ed.ac.uk/software/figtree/).

Figure 1 Mitochondrial genomic organizations of Ranidae frogs. Each tRNA gene is represented by the standard one-letter amino acid code,and S1 = trnSUCN, S2 = trnSAGY, L1 = trnLCUN, L2 = trnLUUR. Other genes are abbreviated as follows: 12S and 16S, 12S and 16S ribosomal RNA;ATP6 and ATP8, adenine triphosphatase subunits 6 and 8; CO1–3, cytochrome c oxidase subunits 1–3; CYTB, cytochrome b; ND1–6 and 4L,NADH dehydrogenase subunits 1–6 and 4L. OL, CR, Ψ, and gap denote replication origin of light strand, D-loop region, pseudogene, and intengenic spacer region, respectively. Genes encoded by the heavy and light strand are denoted at the top and bottom of each gene rectangle box, respectively. The sizes of the boxes do not reflect actual gene length.

3. Results

3.1. Mitogenome Characterization and analysis of three new Rana mitogenomes

3.1.1. Genome organizationThe complete nucleotide sequences of theR. chaochiaoensis,R. kukunorisandR. omeimontismitogenomes have been determined successfully in this study and submitted to the GenBank database under accession numbers KU246048–KU246050 (Table 1). All three mitogenomes were circular, consisting of two rRNAs, 13 PCGs, 22 tRNAs and four intergenic spacer regions (Table S2; Figure S2).The largest intergenic spacer region was located betweenCYTBandtrnLCUN, which was the typical position of D-loop region. We determined the smaller one located in theWANCY-trncluster as OLregion based on its typical stem-loop structure and the surrounding 5’-GCCGG-3’motif (on the light strand). The remaining two gaps were discovered at the two flanks ofND5gene (Figure S2).All three mitogenomes retained the identical genomic organization (Figure 1; Figure S2), and they were 18 591 bp, 18 863 bp, and 19 934 bp in size, respectively (Table 1). The overall base composition of the light strand was 28.85%–29.51% for T, 28.04%–28.45% for C, 27.46%–27.88% for A and 14.56%–15.06% for G with an A + T bias (56.49%–57.39%).

3.1.2. Ribosomal RNA and Protein-Coding genesThe12Sand16S rRNAof three mitogenomes were located betweentrnFandtrnLUURand separated bytrnV. The size of12Sand16S rRNAwere 931 bp and 1582 bp forR.omeimontis, 930 bp and 1576 bp forR. chaochiaoensis,and 929 bp and 1 576 bp forR. kukunoris, respectively(Table S2). The overall base composition of two rRNAs were shown as A ＞ C ＞ T ＞ G.

All mitochondrial genomes shared a set of 13 PCGs,includingND1–6,ND4L,CO1–3,ATP8,ATP6andCYTB, and onlyND6was encoded on the L-strain (Table S2; Figure S2). Most PCGs began with the typical ATG codon, exceptingCO1,ATP6andND4Linitiated with GTG, andND1started at ATC (forR. omeimontis) and GTG (R. chaochiaoensisandR. kukunoris). Six PCGs harbored the traditional complete termination codons TAA (ATP8,ND4LandCYTB), AGA (ND5andND6) and AGG (CO1), whereas the remaining seven PCGs used T(Table S2).

3.1.3. Transfer RNA genesExcluding thetrnSAGYgene,the inferred secondary structures of the other 21 tRNAs of the three mitogenomes conform to the common structural features of mitochondrial tRNAs (Table S2; Figure 1).The base mutations of tRNAs among three mitogenomes existed in the stems and the loops structure.

3.2. Molecular phylogenetic analysisThe final concatenated mtDNA sequence matrix for 48 species was 13 737 bp in size, including 8 777 variable sites of which 974 were singleton sites. Two phylogenetic reconstruction methods (ML and BI) yielded identical tree topologies based on 13 PCGs and two rRNAs, and they favored the following clades and/or relationships of Ranidae (Figure 2): (1) the most basal position of the genusAmnirana;(2) the secondary basal position of the genusGlandirana;(3) the clade ofPelophylax+Amolops; (4) the paraphyly ofBabinainterweaved withSylvirana; (5) the clade ofOdorrana; (6) the monophyly ofLithobatesandRana;(7) the cladeBabina+ (Odorrana+ (Lithobates+Rana)).Within the lineage Ranidae, clade 7 formed the sister taxon to clade 3, but no sufficient statistical support existed for this relationship (BP = 41, BPP = 0.90).

3.3. Ranidae gene rearrangement analysisAccording to our comparison of genome organization, we summarized 10 different gene arrangements (Figure 1;Table 2). All rearrangements occurred in both theND4–trnTand thetrnW–CO1regions (Figure 1; Figure 3).

Our results showed that Type B (also termed as the typical Neobatrachian-type arrangement) was the most common type in ranid (or neobatrachian) mitogenomes.AllPelophylaxfrogs and another twoAmolopsfrogs,A.rickettiandA. wuyiensis(namely theA. rickettispecies group), expressed the Type B. Additionally, Type B was the most basic type, and another nine novel types (from Type C to Type K) were derived from it via diverse rearrangement pathways.

Type C was only discovered inAmnirana albolabris(Figure 2; Table 2). In this type, the positions betweentrnAandtrnN-OL-trnCwere exchanged accompanied with the insertions of some non-coding regions and finally yielding the noveltrnW-gap-trnN-OL-trnC-gaptrnA-gap-trnYorder (Figure 3). Type D was unique to theGlandiranafrogs, which was characterized by thetrnSAGYpseudogene next totrnH(Figure 2; Table 2). Type E was shared byAmolops mantzorumspecies group, which was different from the Type B possessed by theA. rickettispecies group in terms of location of the OLstructure(Figure 2; Table 2). The OLwas translocated from the downstream to the upstream position of thetrnA-trnN,and then several non-coding regions were inserted into this block, yielding the distinctivetrnW-gap-OL-gap-trnA-trnN-gap-trnC-trnYorder (Figure 3).

Figure 2 The ML and BI phylogeny trees derived from the concatenated sequences of 13 protein coding genes and two rRNA genes among Ranidae. Numbers above the lines or beside the nodes are rapid bootstrap proportions calculated with 1 000 replicates and Bayesian posterior probabilities, respectively. The different color represents the different genomic rearrangement features of each species.

(Continued Table 1)

Table 2 Frequency of each mitochondrial genome rearrangement type in family Ranidae.

Table 3 The mitochondrial genome types in the nine genera of the family Ranidae.

Type F was the most common type (32.50%) in ranid mitogenomes so far (Table 3). Type F appeared in most ofRana(including our three species), allLithobates, severalBabinaand oneSylviranafrogs(Figure 2; Table 2). Type G was shared by the twoBabinafrogs (Figure 2; Table 2). A large number of gene rearrangements were found in Type G. The variation of gene rearrangement inOdorranawas quite large, and fourOdorranaspecies held the three types (H, I, and J). In all threeOdorranarearrangement types, thetrnHwas translocated to D-loop downstream,forming aHLTPF-trncluster. Moreover, the position exchange betweentrnNand OLwas only discovered in Type I (O. ishikawae) and Type J (O. schmackeri).In particular, the OLregion was triplicate in Type I (O.ishikawae).R. kunyuensisandR. coreanashared the identical arrangement order Type K. Compared with Type F, this type showed more complex variations: one additional D-loop region was inserted into the upstream ofTPF-trncluster, and theND5was translocated from the typicaltrnSAGYdownstream to thetrnLCUNdownstream(Figure 3).

4. Discussion

4.1. Characteristics analysis of the Rana mitogenomesThreeRanamitogenomes shared the identical genomic organization with those ofR.cf.chensinensis,R.dybowskii,R. huanrensisandR. draytonii(Donget al.,2016; Li et al., 2016a; Figure S1), and this genomic organization was similar to the typical Neobatrachiantype (Kurabayashiet al., 2010; Sumidaet al., 2001).The variation of molecular size and base composition of entire genome among all publishedRanamitogenomes were primarily due to the duplication of D-loop region and the variable numbers of tandem repeat element in D-loop region (Donget al., 2016; Li et al., 2016a, b). The incomplete termination codon T frequently appeared in seven PCGs, and it was completed by post-transcriptional polyadenylation (Ojalaet al., 1981).

Figure 3 Putative mechanism of gene rearrangement of the mitochondrial genome according to the duplication and random loss model. The information of each gene or region is the same as those in Figure 1. The solid arrows represent duplication events and the dashed arrows represent random loss events. The green and blue boxes represent duplication regions; the gray and black boxes represent partial loss and complete deletion, respectively.

4.2. Molecular phylogenetic analysisOverall, the genus level phylogeny reconstructed in our study was congruent with the hypotheses from Liet al.(2014) and Buet al.(2016) but conflicted with other results from some researchers (e.g. Cheet al., 2007; Kurabayashiet al.,2010; Niet al., 2016; Pyron and Wiens, 2011; Wienset al., 2009; Xiaet al., 2014). Our trees placedGlandiranaat a more basal position with strong support(BP = 84, BPP = 1.00), which was in agreement with the result of Buet al.(2016) but different from other reports that located theGlandiranain a nested position within the Ranidae phylogenetic tree with weak statistical support(e.g. Cheet al., 2007; Kurabayashiet al., 2010; Niet al., 2016; Xiaet al., 2014). By reviewing the previous work, we found [Babina] and [Lithobates+Rana] had been considered as the sister group ofOdorrana. Using the single gene or very few genes (e.g. two rRNAs), Cheet al.(2007), Kurabayashiet al.(2010), Wienset al.(2009) and Xiaet al.(2014) found [Babina] was the sister group ofOdorrana. Kakehashiet al.(2013) reconstructed the same phylogenetic relationship using two rRNAs and 13 PCGs and proposed plausible explanation according to the probable gene rearrangement mechanisms (see below).

However, our results robustly supported that[Lithobates+Rana] was the sister group ofOdorrana,which was compatible with other studies based on 13 PCGs (e.g. Bu et al., 2016; Liet al., 2014; Niet al., 2016;Xueet al., 2016). Kakehashiet al.(2013) also noted that the genusBabinaspecies formed a monophyletic group (BP = 100). However, theS. guentheriwas nested inBabinaclade in our phylogenetic trees, as previously reported by Niet al.(2016). The taxonomic history ofS.guentheriwas somewhat complicated (Wuet al., 2016).From 1882 to 2010, this species was successively placed into several genera, such asRana,Hylorana,Hylarana,andBoulengerana(see Frost, 2017). Most recently, it has been classified intoSylviranabased on two mitochondrial and four nuclear gene data (Oliveret al., 2015).Nonetheless, more convincing evidence is indispensable for determining the taxonomic status of this frog.

4.3. Extensive gene rearrangement in RanidaeKurabayashiet al.(2010) stated that the diversity of the mitochondrial genome reorganization in ranids was unexpected. In this study, we summarized 10 different gene orders, and found that all rearrangements occurred at theND4–trnFregion and thetrnW–CO1region.In Caudata mitogenomes, the gene rearrangements also appeared at the two regions (Xiaet al., 2010). In Gymnophiona mitogenomes, the gene rearrangements occurred more at thetrnW–CO1region (San Mauroet al., 2006). Liet al.(2010) indicated that the Anura mitogenome rearrangements mainly occurred at the flanks of D-loopregion, the margin of OLstructure and theIQM-trngenes cluster. Moreover, we found many rearranged patterns, such asWAOLOLOLNCY,WOLANCY,WNOLCAYandWAOLNCY, are discovered in some Ranidae mitochondrial genomes (Figure 1; Figure 2). Therefore,we speculated thetrnW–CO1region and theND4–trnFregion should be the hotspots of Ranidae mitochondrial genomes rearrangement.

AllAmolopsmitogenomes analyzed in this study were classified as Type B and Type E, and they were different from the previously determinedA. larutensisrearrangement type (Kurabayashiet al., 2010; Figure 1),implying that theAmolopsgene rearrangements were various. In particular, the OLregion was triplicate inO.ishikawaemitogenome (Type I). The triploidization of the OLwas unique to this frog in Ranidae, but it was also discovered in the mitogenome ofCallulina kreffti(Brevicipitidae), another Neobatrachia frog (Zhanget al.,2013). In addition, the diploidization of the OLwas found in theA. larutensismitogenome (Kurabayashiet al.,2010).

Interestingly,R. kunyuensisandR. coreanashared one additional D-loop region and duplicate D-loop regions was not unique to these two ranids (Liet al., 2016b),because it was also discovered in another Ranidae speciesA. larutensis(Kurabayashiet al., 2010), and other Neobatrachia taxa, such as Afrobatrachia frogs(Kurabayashi and Sumida, 2013), Mantellidae frogs(Kurabayashiet al., 2006, 2008),Rhacophorus schlegelii(Sanoet al., 2005), andHoplobatrachusspp. (Alamet al., 2010; Yuet al., 2012b). Wanget al.(2015) found that the duplicated D-loop regions within one individual were almost identical in the bushtits mitochondrial genomes, and further supposed that homologous recombination occurred between paralogous D-loop regions from different mtDNA molecule was proposed as the most suitable mechanism for concerted evolution of the duplicated D-loop regions. Unfortunately, in this study we cannot speculate the mechanism for thisRanaduplicated D-loop regions.

4.4. Mechanisms and systematic implication of mitochondrial gene rearrangementGenerally, the vertebrate mitochondrial gene rearrangement was relatively rare and random (Xiaet al., 2014). As stated by many scholars, all observed gene rearrangement events of vertebrate mitogenomes could be classified as translocation, inversion, shuffling, deletion, or duplication(Dowtonet al., 2003; Maceyet al., 1997), and gene shuffling was the prevailing gene rearrangement type(Maceyet al., 1997). In our study, only gene translocation and duplication were discovered in these Ranidae mitogenomes, and gene shuffling was more common than gene duplication (Figure 1; Figure 3). Unlike the D-loop region and OLstructure, which tend to gene duplication,the tRNAs genes and PCGs tend to gene shuffling(Figure 1; Figure 3). For the formation of rearrangement types, several different rearrangement mechanisms were proposed, such as the tandem duplication and random loss model (Maceyet al., 1997; Moritz and Brown, 1987), the tandem duplication and non-random loss model (Lavrovet al., 2002), and the intramitochondrial recombination model (Poultonet al., 1993).

Currently, the duplication and random loss model can be used to explain for most of the animal mitogenome reorganization (e.g. Kakehashiet al., 2013; Kurabayashiet al., 2008). In this model, initially, a duplication including a part of the entire genome happened accidentally because of replication errors (either slippedstrand mispairing or inaccurate termination); then,one of the duplicates of the included genes (or noncoding region) was converted into a pseudogene and subsequently excised from the genome through an accumulation of natural mutations (Dowtonet al., 2003;Maceyet al., 1997; Moritz and Brown, 1987). In the present study, the duplication and random loss model also could explain all rearrangement events discovered in the Ranidae mitogenomes (Figure 3), although some of our views were not compatible with previous views(e.g. Kakehashiet al., 2013; Kurabayashiet al., 2010).Additionally, it was almost impossible that the same gene order was generated independently through different pathways among different taxa.

The vertebrate mitochondrial rearrangement was regarded as unique, random, and a generally rare event(Boore and Brown, 1998; Liu and Huang, 2010; Xiaet al., 2010, 2014), and the occurrence of identical gene rearrangements in two or more lineages indicated that this gene rearrangement type was a synapomorphic type and these lineages were derived from a common ancestor (Maceyet al., 1997), although a few convergent or parallel gene rearrangements have been observed in the vertebrate mtDNAs (e.g. Morrisonet al., 2002;Weiet al., 2014). The remarkable mitochondrial gene rearrangement contributes to our understanding of phylogenetic relationships and is now considered as a valuable molecular marker (Boore and Brown, 1998;Kurabayashiet al., 2008, 2010; Maceyet al., 1997), being widely applied to explore the phylogenetic relationships among various groups at different taxonomic levels (e.g.Kakehashiet al., 2013; Kurabayashiet al., 2006, 2010;Liuet al., 2016; San Mauroet al., 2004, 2014; Weiet al., 2014; Xiaet al., 2010; Xueet al., 2016; Zhanget al., 2008, 2009, 2013).

As mentioned above, the previous studies considered[Babina] or [Lithobates+Rana] as the sister taxon ofOdorrana(e.g. Kakehashiet al., 2013; Kurabayashiet al., 2010; Niet al., 2016; Xueet al., 2016). Additionally,Kakehashiet al.(2013) further pointed out that the[Babina+Odorrana] clade shared a common ancestral gene arrangement type.

Alternatively, we proposed another explanation: all taxa, incluingBabina,Sylvirana,Odorrana,Lithobates,andRana, shared a common ancestral gene order Type F (Figure 2; Figure 3), but this order was completely different from the pattern (ND4-trnH-trnS2-ND5-ND6-trnE-CYTB-D-loop-trnH-trnS2-ND5-ND6-trnE-trnL1-trnT-trnP-trnF-12S-trnV-12S) inferred by Kakehashiet al.(2013). Several lineages possessed their distinctive gene rearrangements, includingGlandiranaspp.,Amolops mantzorumspecies group,Amolops rickettispecies group,Pelophylaxspp., and theRana+Lithobateslineage (excludingR. kunyuensisandR. coreana). The genusAmolopswas a complicated group. In siblingA.larutensis, a lot of mitochondrial gene rearrangements(including duplication of D-loop region, duplication of OL, transpositions oftrnK,trnHandtrnG-ND3block)had been discovered by Kurabayashiet al.(2010).Considering the fact that this species possessed a nested position withinAmolops, Kurabayashiet al.(2010)inferred the genomic reorganization was likely to have occurred in a common ancestor ofAmolops, or during the diversification of this taxon. Now, the latter was confirmed by more available mitogenomes. InGlandiranaspp., thetrnSAGYpseudogene was proved as a valuable molecular marker for its phylogenetic analysis.

5. Conclusion

The threeRanafrogs shared the identical mitogenome arrangement type, which was extremely similar to the typical Neobatrachian-type arrangement shared by most frogs. The phylogenetic analysis using PCGs and rRNAs sequences from 55 mitogenomes demonstrated that the genusAmniranaoccupied the most basal position among the Ranidae and the [Lithobates+Rana] was the closest sister group ofOdorrana. The diversity of Ranidae mitogenome arrangements was unexpected high, and the 47 mitogenomes of 40 ranids were classified into 10 different gene rearrangement types. Some taxa owned their distinctive gene rearrangement characteristics,which had significant implication for their phylogeny analysis. The tandem duplication and random loss model can explain all rearrangement events discovered in all Ranidae mitogenomes.

AcknowledgementsWe wish to offer our sincere thanks to everyone who helped make this paper a reality.This work was supported by the Innovative Research Team in University of Sichuan Bureau of Education(No.14TD0002) and the Scientific Research Fund of Sichuan Provincial Education Department (No.11ZA077).

Alam M. S., Kurabayashi A., Hayashi Y., Sano N., Khan M. R.,Fujii T., Sumida M. 2010. Complete mitochondrial genomes and novel gene rearrangements in two dicroglossid frogs,Hoplobatrachus tigerinusandEuphlyctis hexadactylus, from Bangladesh. Genes Genet Syst, 85: 219–232

Anderson S., Bankier A. T., Barrell B. G., de Bruijn M. H. L.,Coulson A. R., Drouin J., Eperon I. C., Nierlich D. P., Roe B.A., Sanger F., Schreier P. H., Smith A. J. H., Staden R., Young I. G. 1981. Sequence and organization of the mitochondrial human genome. Nature, 290: 457–465

AmphibiaWeb. 2017. Information on amphibian biology and conservation. Retrieved from http://amphibiaweb.org

Boore J. L., Brown W. M. 1998. Big trees from little genomes:mitochondrial gene order as a phylogenetic tool. Curr Opin Genet Dev, 8: 668–674

Bu X., Zhang L., He K., Jiang Y., Nie L. 2016. The complete mitochondrial genome of theOdorrana schmackeri(Anura,Ranidae). Mitochondrial DNA Part B, 1: 162–163

Cameron S. L. 2014. How to sequence and annotate insect mitochondrial genomes for systematic and comparative genomics research. Syst Entomol, 39: 400–411

Che J., Pang J., Zhao H., Wu G. F., Zhao E. M., Zhang Y. P.2007. Phylogeny of Raninae (Anura: Ranidae) inferred from mitochondrial and nuclear sequences. Mol Phylogenet Evol, 43:1–13

Chen G., Wang B., Liu J., Xie F., Jiang J. 2011. Complete mitochondrial genome ofNanorana pleskei(Amphibia:Anura: Dicroglossidae) and evolutionary characteristics of the amphibian mitochondrial genomes. Curr Zool, 57: 785–805

Chen Z., Zhang J., Zhai X., Zhu Y., Chen X. 2015. Complete mitochondrial genome of the green odorous frogOdorrana margaretae(Anura: Ranidae). Mitochondrial DNA Part A, 26:487–488

Darriba D., Taboada G. L., Doallo R., Posada D. 2012.jModelTest 2: More models, new heuristics and parallel computing. Nat Methods, 9: 772

Desjardins P., Morais R. 1990. Sequence and gene organization of the chicken mitochondrial genome: A novel gene order in higher vertebrates. J Mol Biol, 212: 599–634

Dong B., Yu Z., Yang B. 2016. The complete mitochondrial genome of theRana huanrensis(Anura: Ranidae). Mitochondrial DNA Part A, 27: 4551–4552

Dowton M., Castro L. R., Campbell S. L., Bargon S. D., Austin A. D. 2003. Frequent mitochondrial gene rearrangements at the hymenopteran nad3–nad5 junction. J Mol Evol,56: 517–526

Frost D. R. 2017. Amphibian Species of the World: an Online Reference Version 6.0. Retrieved from http://research.amnh.org/herpetology/amphibia/index.html

Guindon S., Dufayard J. F., Lefort V., Anisimova M., Hordijk W., Gascuel O. 2010. New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0. Syst Biol, 59: 307–321

Kakehashi R., Kurabayashi A., Oumi S., Katsuren S., Hoso M.,Sumida M. 2013. Mitochondrial genomes of Japanese Babina frogs (Ranidae, Anura): unique gene arrangements and the phylogenetic position of genusBabina. Genes Genet Syst, 88:59–67

Kumazawa Y., Endo H. 2004. Mitochondrial genome of the Komodo dragon: efficient sequencing method with reptileoriented primers and novel generearrangements. DNA Res, 11:115-125

Kurabayashi A., Sumida M. 2013. Afrobatrachian mitochondrial genomes: Genome reorganization, gene rearrangement mechanisms, and evolutionary trends of duplicated and rearranged genes. BMC Genomics, 14: 633

Kurabayashi A., Sumida M., Yonekawa H., Glaw F., Vences M.,Hasegawa M. 2008. Phylogeny, recombination, and mechanisms of stepwise mitochondrial genome reorganization in mantellid frogs from Madagascar. Mol Biol Evol, 25: 874–891

Kurabayashi A., Usuki C., Mikami N., Fujii T., Yonekawa H.,Sumida M., Hasegawa M. 2006. Complete nucleotide sequence of the mitochondrial genome of a Malagasy poison frogMantella madagascariensis: evolutionary implications on mitochondrial genomes of higher anuran groups. Mol Phylogenet Evol, 39:223–236

Kurabayashi A., Yoshikawa N., Sato N., Hayashi Y., Oumi S., Fujii T., Sumida M. 2010. Complete mitochondrial DNA sequence of the endangered frogOdorrana ishikawae(family Ranidae) and unexpected diversity of mt gene arrangements in ranids. Mol Phylogenet Evol, 56: 543–553

Larkin M. A., Blackshields G., Brown N. P., Chenna R.,McGettigan P. A., McWilliam H., Valentin F., Wallace I. M.,Wilm A., Lopez R., Thompson J. D., Gibson T. J., Higgins D.G. 2007. Clustal W and Clustal X version 2.0. Bioinformatics,23: 2947–2948

Lavrov D. V., Boore J. L., Brown W. M. 2002. Complete mtDNA sequences of two millipedes suggest a new model for mitochondrial gene rearrangements: duplication and nonrandom loss. Mol Biol Evol, 19: 163–169

Li J., Lei G., Fu C. 2016a. Complete mitochondrial genomes of two brown frogs,Rana dybowskiiandRanacf.chensinensis(Anura: Ranidae). Mitochondrial DNA Part A, 27: 155–156

Li J., Yin W., Xia R., Lei G., Fu C. 2016b. Complete mitochondrial genome of a brown frog,Rana kunyuensis(Anura: Ranidae).Mitochondrial DNA Part A, 27: 34–35

Li X. Q., Zhang M., Wang Y. S., Wu X. B. 2010. Research rogress on mitochondrial genomes and gene rearrangements in Anura. In: Ji X. (Ed) Herpetological sinaca 12thedn. Southeast University Press, Nanjing, pp: 387–394

Li Y., Zhang H., Wu X., Xue H., Yan P., Wu X. 2014. A novel mitogenomic rearrangement forOdorrana schmackeri(Anura:Ranidae) and phylogeny of Ranidae inferred from thirteen mitochondrial protein-coding genes. Amphibia-Reptilia, 35:331–343

Liu J., Xue R., Wang Y., Li D., Yan Q., Yang J. 2016. The nearcomplete mitogenome sequence of the Omei Horned ToadMegophrys omeimontisLiu, 1950 (Anura, Megophryidae).Mitochondrial DNA Part A, 27: 2389–2390

Liu N., Huang Y. 2010. Complete mitochondrial genome sequence ofAcrida Cinerea(Acrididae: Orthoptera) and comparative analysis of mitochondrial genomes in orthoptera. Int J Genomics,2010: 319486

Liu Z. Q., Wang Y. Q., Su B. 2005. The mitochondrial genome organization of the rice frog,Fejervarya limnocharis(Amphibia:Anura): a new gene order in the vertebrate mtDNA. Gene, 346:145–151

Luo A., Qiao H., Zhang Y., Shi W., Ho S. Y. W., Xu W.,Zhang A., Zhu C. 2010. Performance of criteria for selecting evolutionary models in phylogenetics: a comprehensive study based on simulated datasets. BMC Evol Biol, 10: 242

Macey J. R., Larson A., Ananjeva N. B., Fang Z., Papenfuss T. J. 1997. Two novel gene orders and the role of light-strand replication in rearrangement of the vertebrate mitochondrial genome. Mol Biol Evol, 14: 91–104

Mindell D. P., Sorenson M. D., Dimcheffe D. E. 1998. Multiple independent origins of mitochondrial gene order in birds. Proc Natl Acad Sci USA, 95: 10693–10697

Moritz C., Brown W. M. 1987. Tandem duplications in animal mitochondrial DNAs: variation in incidence and gene content among lizards. Proc Natl Acad Sci USA, 84: 7183–7187

Morrison C. L., Harvey A. W., Lavery S., Tieu K., Huang Y.,Cunningham W. C. 2002. Mitochondrial gene rearrangements confirm the parallel evolution of the crab-like form. Proc R Soc Lond B, 269: 345–350

Mueller R. L., Boore J. L. 2005. Molecular mechanisms of extensive mitochondrial gene rearrangement in plethodontid salamanders. Mol Biol Evol, 22: 2104–2112

Ni N., Yu D., Storey K. B., Zheng R., Zhang J. 2016. The complete mitochondrial genome ofLithobates sylvaticus(Anura:Ranidae). Mitochondrial DNA Part A, 27: 2460–2461

Ojala D., Montoya J., Attardi G. 1981. tRNA punctuation model of RNA processing in human mitochondria. Nature, 290: 470–474

Oliver L. A., Prendini E., Kraus F., Raxworthy C. J. 2015.Systematics and biogeography of theHylaranafrog (Anura:Ranidae) radiation across tropical Australasia, Southeast Asia,and Africa. Mol Phylogenet Evol, 90: 176–192

Pabijan M., Spolsky C., Uzzell T., Szymura J. M. 2008.Comparative analysis of mitochondrial genomes inBombina(Anura; Bombinatoridae). J Mol Evol, 67: 246–256

Poulton J., Deadman M. E., Bindoff L., Morten K., Land J.,Brown G. 1993. Families of mtDNA re-arrangements can be detected in patients with mtDNA deletions: duplications may be a transient intermediate form. Hum Mol Genet, 2: 23–30

Pyron R. A., Wiens J. J. 2011. A large-scale phylogeny of Amphibia including over 2800 species, and a revised classification of extant frogs, salamanders, and caecilians. Mol Phylogenet Evol, 61: 543–583

Roe B. A., Ma D. P., Wilson R. K. E., Wong J. F. H. 1985.The complete nucleotide sequence of theXenopus laevismitochondrial genome. J Biol Chem, 260: 9759–9774

Ronquist F., Huelsenbeck J. P. 2003. MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics, 19:1572–1574

San Mauro D., Gower D. J., Oommen O. V., Wilkinson M., Zardoya R. 2004. Phylogeny of caecilian amphibians(Gymnophiona) based on complete mitochondrial genomes and nuclear RAG1. Mol Phylogenet Evol, 33: 413–427

San Mauro D., Gower D. J., Zardoya R., Wilkinson M. 2006. A hotspot of gene order rearrangement by tandem duplication and random loss in the vertebrate mitochondrial genome. Mol Biol Evol, 23: 227–234

San Mauro D., Gower D. J., Müller H., Loader S. P., Zardoya R., Nussbaum R. A., Wilkinson M. 2014. Life-history evolution and mitogenomic phylogeny of caecilian amphibians.Mol Phylogenet Evol, 73: 177–189

Sano N., Kurabayashi A., Fujii T., Yonekawa H., Sumida M.2005. Complete nucleotide sequence of the mitochondrial genome of Schlegel's tree frogRhacophorus schlegelii(family Rhacophoridae): duplicated control regions and gene rearrangements. Genes Genet Syst, 80: 213–224

Su X., Wu X. B., Yan P., Cao S. Y., Hu Y. L. 2007. Rearrangement of a mitochondrial tRNA gene of the concave-eared torrent frog,Amolops tormotus. Gene, 394: 25–34

Sumida M., Kanamori Y., Kaneda H., Kato Y., Nishioka M.,Hasegawa M., Yonekawa H. 2001. Complete nucleotide sequence and gene rearrangement of the mitochondrial genome of the Japanese pond frogRana nigromaculata. Genes Genet Syst, 76: 311–325

Sumida M., KondoY., Kanamori Y., Nishioka M. 2002. Interand intraspecific evolutionary relationships of the rice frogRana limnocharisand the allied speciesR. cancrivorainferred from crossing experiments and mitochondrial DNA sequences of the 12S and 16S rRNA genes. Mol Phylogenet Evol, 25: 293–305

Tamura K., Stecher G., Peterson D., Filipski A., Kumar S. 2013.MEGA6: molecular evolutionary genetics analysis version 6.0.Mol Biol Evol, 30: 2725–2729

Tzeng C. S., Hui C. F., Shen S. C., Huang P. C. 1992. The complete nucleotide sequence of theCrossostoma lacustremitochondrial genome: conservation and variations among vertebrates. Nucl Acids Res, 20: 4853–4858

Wang X., Huang Y., Liu N., Yang J., Lei F. 2015. Seven complete mitochondrial genome sequences of bushtits (Passeriformes,Aegithalidae,Aegithalos): the evolution pattern in duplicated control regions. Mitochondrial DNA, 26: 1–7

Wei S. J., Li Q., Achterberg K. V., Chen X. X. 2014. Two mitochondrial genomes from the families Bethylidae and Mutillidae: independent rearrangement of protein-coding genes and higher-level phylogeny of the Hymenoptera. Mol Phylogenet Evol, 77: 1–10

Wiens J. J., Sukumaran J., Pyron R. A., Brown R. M. 2009.Evolutionary and biogeographic origins of high tropical diversity in Old World frogs (Ranidae). Evolution, 63: 1217–1231

Wu X., Li Y., Zhang H., Jiang Z., Xue H., Yan P., Wu X. 2016.The complete mitochondrial genome ofHylarana guentheri(Amphidia, Anura, Ranidae). Mitochondrial DNA Part A, 27:1223–1224

Xia Y., Peng R., Zeng X. M. 2010. Characteristics of mitochondrial gene rearrangement in Caudata. In: Ji X. (ed) Herpetological sinaca 12thedn. Southeast University Press, Nanjing, pp: 363–371

Xia Y., Zheng Y., Miura I., Wong P. B., Murphy R. W., Zeng X.2014. The evolution of mitochondrial genomes in modern frogs(Neobatrachia): nonadaptive evolution of mitochondrial genome reorganization. BMC Genomics, 15: 691

Xue R., Liu J., Yu J., Yang J. 2016. The complete mitogenome ofAmolops loloensisand related phylogenetic relationship among Ranidae. Mitochondrial DNA Part A, 27: 4629–4630

Yan L., Geng Z. Z., Yan P., Wu X. B. 2016. The complete mitochondrial genome ofGlandirana tientaiensis(Ranidae,Anura). Mitochondrial DNA Part A, 27: 1154–1155

Yu D., Zhang J., Zheng R. 2012a. The complete mitochondrial genome ofBabina adenopleura(Anura: Ranidae). Mitochondrial DNA Part A, 23: 423–425

Yu D., Zhang J., Zheng R., Shao C. 2012b. The complete mitochondrial genome ofHoplobatrachus rugulosus(Anura:Dicroglossidae). Mitochondrial DNA Part A, 23: 336–337

Zardoya R., Garrido-Pertierra A., Bautista J. M. 1995. The complete nucleotide sequence of the mitochondrial DNA genome of the rainbow trout,Oncorhynchus mykiss. J Mol Evol,41: 942–951

Zhang P., Papenfuss T. J., Wake M. H., Qu L., Wake D. B.2008. Phylogeny and biogeography of the family Salamandridae(Amphibia: Caudata) inferred from complete mitochondrial genomes. Mol Phylogenet Evol, 49: 586–597

Zhang P., Wake M. H. 2009. A mitogenomic perspective on the phylogeny and biogeography of living caecilians (Amphibia:Gymnophiona). Mol Phylogenet Evol, 53: 479–491

Zhang P., Liang D., Mao R. L., Hillis D. M., Wake D. B.,Cannatella D. C. 2013. Efficient sequencing of anuran mtDNAs and a mitogenomic exploration of the phylogeny and evolution of frogs. Mol Biol Evol, 30: 1899–1915