Complete mitochondrial genomes of Nanorana taihangnica and N. yunnanensis (Anura: Dicroglossidae) with novel gene arrangements and phylogenetic relationship of Dicroglossidae

Background Complete mitochondrial (mt) genomes have been used extensively to test hypotheses about microevolution and to study population structure, phylogeography, and phylogenetic relationships of Anura at various taxonomic levels. Large-scale mt genomic reorganizations have been observed among many fork-tongued frogs (family Dicroglossidae). The relationships among Dicroglossidae and validation of the genus Feirana are still problematic. Hence, we sequenced the complete mt genomes of Nanorana taihangnica (=F. taihangnica) and N. yunnanensis as well as partial mt genomes of six Quasipaa species (dicroglossid taxa), two Odorrana and two Amolops species (Ranidae), and one Rhacophorus species (Rhacophoridae) in order to identify unknown mt gene rearrangements, to investigate the validity of the genus Feirana, and to test the phylogenetic relationship of Dicroglossidae. Results In the mt genome of N. taihangnica two trnM genes, two trnP genes and two control regions were found. In addition, the trnA, trnN, trnC, and trnQ genes were translocated from their typical positions. In the mt genome of N. yunnanensis, three control regions were found and eight genes (ND6, trnP, trnQ, trnA, trnN, trnC, trnY and trnS genes) in the L-stand were translocated from their typical position and grouped together. We also found intraspecific rearrangement of the mitochondrial genomes in N. taihangnica and Quasipaa boulengeri. In phylogenetic trees, the genus Feirana nested deeply within the clade of genus Nanorana, indicating that the genus Feirana may be a synonym to Nanorana. Ranidae as a sister clade to Dicroglossidae and the clade of (Ranidae + Dicroglossidae) as a sister clade to (Mantellidae + Rhacophoridae) were well supported in BI analysis but low bootstrap in ML analysis. Conclusions We found that the gene arrangements of N. taihangnica and N. yunnanensis differed from other published dicroglossid mt genomes. The gene arrangements in N. taihangnica and N. yunnanensis could be explained by the Tandem Duplication and Random Loss (TDRL) and the Dimer-Mitogenome and Non-Random Loss (DMNR) models, respectively. The invalidation of the genus Feirana is supported in this study.

In the present study, we determined the complete mt genomes of N. taihangnica and N. yunnanensis as well as the partial mt genomes of six Quasipaa species (Dicroglossidae), two Odorrana and two Amolops species (Ranidae), and one Rhacophorus species (Rhacophoridae). In this paper, we follow the system of anuran taxonomy published by Fei et al. [42] and Frost et al. [43] to prevent unnecessary confusion in taxonomy. The data was used to determine unknown mt gene rearrangements, to investigate the validity of the genus Feirana, and to test the phylogenetic relationships of Ranidae and Dicroglossidae.

Ethical statement
The thirteen species studied (N. taihangnica, N. yunnanensis, Q. boulengeri, Q. exilispinosa, Q. jiulongensis, Q. robertingeri, Q. shini, Q. verrucospinosa, Odorrana livida, O. schmackeri, Amolops hongkongensis, A. wuyiensis, Rhacophorus dennysi) are not protected by the provisions of the laws of People's Republic of China on the protection of wildlife. Thus, the experiments in this study were performed with toe-clip tissue samples collected from all frog specimens and stored in 100% ethanol. Sample acquisition was reviewed, approved and carried out in accordance with the relevant guidelines of the Committee of Animal Research Ethics of Zhejiang Normal University.

Sample collection
Specimens included two species of Nanorana (N. taihangnica = F. taihangnica, N. yunnanensis), seven samples belonging to six species of Quasipaa (Q. boulengeri, Q. exilispinosa, Q. jiulongensis, Q. robertingeri, Q. shini, Q. verrucospinosa) (Dicroglossidae) including two Q. boulengeri samples from two different sites, two species of Odorrana (O. livida, O. schmackeri) (Ranidae), two species of Amolops (A. hongkongensis, A. wuyiensis) (Ranidae), and R. dennysi (Rhacophoridae). Information on all of the sequenced samples is shown in Table 1. We were unable to successfully sequence the displacement loop (D-loop) region of these samples except for N. taihangnica and N. yunnanensis because of highly repetitive regions in the D-loop or other unknown reasons despite many optimization efforts; this is similar to the report of Zhang et al. [45].

PCR and sequencing
Total DNA was extracted from the clipped toe of each frog specimen using a DNeasy Tissue Kit (Qiagen, Germany). We amplified overlapping fragments that covered the entire mt genome of N. taihangnica and N. yunnanensis by normal PCR and long-and-accurate polymerase chain reaction (LA PCR) methods slightly modified from Yu et al. [5,46] and Zhang et al. [45]. All PCR procedures were performed using a MyCycler Thermal Cycler (Bio-Rad, Hercules, CA, USA). TaKaRa Ex-Taq and LA-Taq kits (Takara Biomedical, Dalian, China) were used for the normal and LA-PCR reactions. The resulting PCR fragments were electrophoresed on 1% agarose gels, and all target DNAs were purified from excised pieces of gel using a SanPrep DNA Gel Extraction Kit (Sangon Biotech, Shanghai, China) prior to sequencing. The sequences for each fragment were obtained in an automated DNA sequencer (ABI 3730) from both strands. The long fragments were sequenced using specific primer walking of both strands.

Sequence assembly and analysis
Sequences were checked and assembled using SeqMan (Lasergene version 5.0) [47]. The locations of the 13 protein coding genes and two rRNA genes were determined by comparison with the available RefSeq sequences of closely related anurans downloaded from GenBank using ClustalW in Mega 5.0 [48,49]. All tRNA genes were identified by their cloverleaf secondary structure using tRNAscan SE 1.21 [50] or determined by comparison with the homologous sequences of other anurans. The mt genomes (see Fig. 1) of all taxa were analyzed to determine the corresponding mt gene arrangements. The resultant sequences were deposited in GenBank with accession numbers KF199146-KF199152, KX233864-KX233869 and KM282625 (see Table 1).

Molecular phylogenetic analysis
With the recently increased number of mitochondrial genomes available for Anura, phylogenetic analyses were performed with 83 anurans for which complete or partial mt genomes were available including 14 samples of the 13 species from this study. In total this included the ingroup of 33 species from Ranidae [1,27,28,45,46,[51][52][53][54][55][56][57][58][59][60][61][62][63], 28 species from Dicroglossidae [5, 10, 12-14, 16, 17, 31-33, 64], 13 species from Rhacophoridae [11,18,65], one species from Mantellidae [9], one species from Petropedetidae [45], one species from Pyxicephalidae [45], one species from Phrynobatrachidae [45], one species from Ptychadenidae [45], one species from Brevicipitidae (outgroup) [26], one species from Hyperoliidae (outgroup) [26] and two species of Microhylidae (outgroups) [66,67]. In order to discuss the phylogenetic relationship of Anura, we used the amino acid data and the nucleotides data to compare the identical topology or not according to the methods of Zhang et al. [6] and Zhou et al. [12]. The amino acid sequences of 10 mt protein-coding genes were separately aligned in Mega 5.0 [48] excluding the ATP8, ND5 and ND6 genes for the following reasons: (a) the ATP8 sequence was too short in length and had too little good information (only 18 nucleotides or < 0.5% of the total nucleotides of combined PCGs) after G-Block analysis, (b) the loss of the ND5 gene in some species [18], and (c) the heterogeneous base composition and poor phylogenetic performance for ND6 which failed to support the consistency analysis with other PCGs [45]. The alignments were revised using Gblocks 0.91b software with the default parameters [68] to select conserved regions of the putative amino acids. We concatenated the alignments of the 10 other mitochondrial protein-coding genes and got an alignment consisting of 2497 amino acid residues as 10Paa dataset. An alignment of 7491 nucleotides sites with 4919 variable informative sites was converted from 2497 amino acids data directly using the amino acid alignment as the backbone. Saturation analysis was performed for subsets with first, second, and third codon positions using DAMBE 4.2.13 [69]. The results showed that the third codon positions were saturated. Thus, we excluded the third codon positions from further phylogenetic analyses and obtained a dataset called 10P consisting of 4994 nucleotide sites from the 1st and 2nd codon positions of the 10 protein-coding genes according to the methods of Cameron et al. [70], Zhang et al. [6] and Zhou et al. [12]. The phylogeny was analyzed using the combined datasets 10P (nucleotides dataset) and 10Paa (amino acid dataset) by the maximum likelihood (ML) and Bayesian inference (BI) methods. To improve the fit of the substitution model to the datasets of 10P and 10Paa, we compared data partitioning schemes according to the Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) using the program PartitionFinder v1.0 and PartitionFinderProtein [71]. We set the 10 coding-  genes as 20 partitions in dataset 10P and 10 partitions in dataset 10Paa, respectively. For the dataset 10P, twentypartitions were optimal: 1) first codon positions of the 10 protein-coding genes; 2) second codon positions of the 10 protein-coding genes. The best substitution model of twenty-partitions in ten different genes of dataset 10P is always GTR + I + G. For the dataset 10Paa, ten-partitions were optimal: 10 protein-coding genes with MTMAM. So, the optimal model for 10P with twenty partitions and the optimal model for 10Paa with ten partitions was chosen for ML by the RaxML program [72] and Bayesian analyses by MrBayes3.1.2 [73][74][75], respectively. ML and BI analyses for datasets of 10P and 10Paa were separately performed using the RaxML program [72] with 1000 bootstrap replications and a modified version of MrBayes3.1.2 [73]. During BI analysis, the following settings were applied: number of Markov chain Monte Carlo (MCMC) generations = 10 million; sampling frequency = 1000; burn-in = 1000. The burn-in size was determined by checking convergences of -log likelihood (−ln L). The robustness of the resulting ML tree was evaluated using bootstrap percentages calculated from nonparametric bootstrap analyses, and statistical support of the resulting BI trees was determined based on Bayesian posterior probability (BPP).

Genome organization of mtDNA
The N. taihangnica mt genome is 21,322 base pairs (bp) in length and contains 13 protein coding genes, two rRNA genes, 24 tRNA genes (including extra trnM and trnP genes), and 10 NCRs including two control regions (CRs). The two CRs were located between the cytochrome b (Cyt b) and trnL genes (CR1 2014 bp) and between the trnC and trnT genes (CR2 2698 bp).
Remarkably, CR1 and CR2 have nearly identical nucleotide sequences (99.9% similarity with only 1 substitution in 2014 alignment sites) excluding the extra 5′-635 bases and 3′-49 bases in CR2. Tandem duplication of the trnM gene and an additional trnP gene were found (Fig. 1).
The trnT-trnP-trnF tRNA cluster moved from the typical neobatrachian LTPF tRNA cluster to a position between the CR1 and NADH dehydrogenase subunit 1 (ND1) genes. The typical LTPF tRNA cluster was replaced by a trnL-trnP-pseudo trnF tRNA cluster. The pseudo-trnF showed 89.9% nucleotide similarity with the corresponding trnF gene in the trnT-trnP-trnF tRNA cluster. This pseudo-trnF contained the same anticodon nucleotides (Fig. 2) compared to trnF. The trnA, trnN, trnC, and trnQ genes were translocated from their typical positions and replaced by a 40-138 bp NCR (Fig. 1). The trnQ gene moved from the typical dicroglossid IQMM tRNA cluster to a location between a 209 bp NCR and a 208 bp NCR (Fig. 1) and within the former IQMM tRNA cluster the trnQ gene was replaced by a 40 bp NCR between the trnI and tandem trnM genes. The trnA, trnN, and trnC genes also moved from the WANCY tRNA cluster to a position between a 208 bp NCR and CR1 (Fig. 1). The positions of the Light-strand replication origin (O L ) are located between a 138 bp NCR (non-coding region) and a 52 bp NCR for the translocations of trnA, trnN, and trnC genes, and a W-NCR (138 bp)-O L -NCR (52 bp)-Y gene cluster was formed in the position of the typical WANO L CY gene cluster. Furthermore, a new cluster consisting of a L-NCR (209 bp)-Q-NCR (208 bp)-A-N-C gene arrangement was observed (Fig. 1).
The two trnP genes contained the same anticodon nucleotides (Fig. 2). The trnF pseudogene contained the same anticodon nucleotides as in trnF whereas trnA and trnN contained different anticodon nucleotides (Fig. 2). The N. yunnanensis mt genome is 23,685 bp in length and contains 13 protein coding genes, two rRNA genes, 23 tRNA genes (including an extra trnM gene), and nine non-coding regions (including three control regions) (Fig. 1). Eight genes (ND6, trnP, trnQ, trnA, trnN, trnC, trnY and trnS genes) in the L-stand were translocated from the typical position to CR regions or near to CR regions and grouped together. CR1, CR2 and CR3 with lengths of 1635 bp, 1581 bp and 1560 bp, respectively, were found between Cyt b and trnQ, between ND6 and trnP, and between trnS and trnL, respectively (Fig. 1). The three CRs have a similar sequence with a length of 1372 bp. The typical WAN(O L )CY tRNA cluster was replaced by a modified W-NCR (139 bp) -O L -NCR (84 bp) arrangement. Through trnP translocation, the LTF tRNA cluster replaced the LTPF tRNA cluster. Through trnQ translocation, the I-NCR (40 bp)-MM tRNA cluster replaced the IQMM tRNA cluster. Through ND6 gene translocation, a 231 bp NCR replaced the ND6 gene in the original region. Through translocation of the trnS gene to between cytochrome c oxidase subunit I (COI) and trnD, the 51 bp NCR replaced trnS. A 47 bp NCR was found between the trnS and ND5 genes.
The detailed gene rearrangements of other known dicroglossids, ranids and rhacophorids in this study are described below.
Quasipaa boulengeri, Q. jiulongensis, Q. verrucospinosa The typical WANO L CY tRNA cluster was replaced by a W-A-NCR-O L -NCR-N-NCR-C-Y tRNA cluster, the NCR of which ranged from 20 bp to 201 bp. The typical IQMM tRNA cluster was found. A 41-56 bp NCR was also found between the trnS and ND5 genes.

Quasipaa robertingeri
The typical WANO L CY tRNA cluster was replaced by a W-A-N-310 bp NCR-O L -C-Y tRNA cluster. The typical IQMM tRNA cluster was found. A 46 bp NCR was also found between the trnS and ND5 genes.

Quasipaa exilispinosa, Q. shini
The typical WANO L CY and IQMM tRNA clusters were also found in Q. spinosa and Q. yei. A 32-48 bp NCR was found between the trnS and ND5 genes.

Odorrana livida
The typical WANO L CY tRNA cluster and a 52 bp NCR between the ND5 and ND6 gene were found.

Odorrana schmackeri
The typical WANO L CY tRNA cluster was replaced by a W-A-NCR-O L-NCR-C-Y tRNA cluster. A 296 bp NCR with tandem sequence between the ND5 and ND6 genes was found.

Amolops hongkongensis, A. wuyiensis
The typical WANO L CY tRNA cluster and no NCR between the ND5 and ND6 gene were found.

Rhacophorus dennysi
The typical WANO L CY tRNA cluster was found and the ND5 gene between trnS and ND6 was translocated to the region between CR and trnT.

Phylogenetic analysis
All BI and ML phylogenetic analyses performed in this study showed similar topologies (Figs. 3 and 4). In the phylogeny of Dicroglossidae, Ranidae, Mantellidae and Rhacophoridae, the monophyly of Dicroglossidae, Ranidae and Rhacophoridae are well supported. Dicroglossidae is a sister clade of Ranidae (1.00 in posterior probability for nucleotides and amino acids datasets; 64% and 69% bootstrap frequencies for nucleotides and amino acids, respectively), and Mantellidae is a sister clade of Rhacophoridae (1.00 in posterior probability for nucleotides and amino acids datasets; 58% and 75% bootstrap frequencies for nucleotides and amino acids, respectively). Then the clade of (Ranidae + Dicroglossidae) is a sister clade of (Mantellidae + Rhacophoridae) (1.00 in posterior probability for nucleotides and amino acids datasets; 46% and 59% bootstrap frequencies for nucleotides and amino acids, respectively).

Discussion
The mtDNA arrangement In Dicroglossidae, a 32-85 bp NCR between the trnS (AGY) and ND5 genes was observed in Paini and    Limnonectini, while a 7-39 bp NCR between the trnS (AGY) and ND6 genes was observed in Dicroglossini and Occidozygini for the translocation of the ND5 gene (Fig. 1). The presence of short non-coding sequences among the rearranged genes has also been observed in previous studies [5]. In this study, both N. yunnanensis and N. taihangnica have a short apomorphic NCR between trnI and trnM; this was also reported in N. quadranus by Zhang et al. [76].
Interspecies tRNA gene rearrangements are well known [11][12][13][14][15][16], but few were found in the current work. Comparing the tRNA gene rearrangements of mt genomes in all known frogs, we found that the LTPF tRNA cluster, the IQMM tRNA cluster and the WANCY tRNA cluster can easily undergo gene rearrangement, the phenomenon appearing not only in interspecies but also intraspecies comparisons (eg. N. taihangnica and Q. boulengeri). Comparing two mt genomes of N. taihangnica (=F. taihangnica) between this study and a previously sequenced N. taihangnica [13], we found that different gene rearrangements of the trnP, trnF, trnQ, trnA, trnN and trnC genes, the IQMM tRNA cluster and the WANCY tRNA cluster existed. The trnT gene of the LTPF tRNA cluster was lost in the previously sequenced N. taihangnica [13], whereas the trnT gene between CR1 and trnP gene was found in N. taihangnica of this study. The L-NCR (35 bp) -PF tRNA cluster in the previously sequenced N. taihangnica [13] was also found in N. taihangnica of this study but an extra TPF tRNA cluster between CR1 and a 289 bp NCR occurred in N. taihangnica of this study. The IQMM tRNA cluster was found in previously sequenced N. taihangnica whereas the I-NCR (40 bp)-MM tRNA cluster was found in N. taihangnica of this study because the trnQ gene was translocated. The WANO L CY tRNA cluster was found in previously sequenced N. taihangnica while W-NCR (138 bp)-O L -NCR (52 bp)-Y was found in N. taihangnica of this study because of the translocation of trnA, trnN and trnC. Comparing mt genomes of Q. boulengeri between this study and other known sequences [16,30], we found different gene rearrangements in the WANO L CY tRNA cluster as also found by Xia et al. [77]. In species of the Genus Nanorana and Quasipaa, two types of tRNA clusters (I-NCR-MM or IQMM, WANO L CY or WAO L NCY) were found. Even in the same species, N. taihangnica and Q. boulengeri, different tRNA clusters were found, which may motivate future discussions on mitochondrial gene arrangements among Nanorana and Quasipaa species. This suggests that more mt genomes of Nanorana and Quasipaa species need to be sequenced to further determine how these different gene arrangements formed.
In N. yunnanensis, seven tRNA genes (trnQ, trnA, trnC, trnY, trnS, trnN and trnP) and the ND6 gene on the L-stand were translocated into or near to control regions and grouped together. We did not find any other species of Anura where these gene arrangements existed, but in a fish Crossorhombus azureus (Pleuronectiformes: Bothidae) [23] seven tRNA genes (trnQ, trnA, trnC, trnY, trnS1, trnE, trnP) and the ND6 gene encoded by the Light-strand (L-strand) were translocated to a position between trnT and trnF, which is very similar to our results.

Possible gene rearrangement mechanisms
In N. taihangnica, we observed several gene rearrangements (extra CR, trnP and trnM genes as well as translocation of trnQ, trnA, trnN, and trnC) in the region between CR and the cox1 gene. We propose that the gene rearrangements may be explained by the TDRL model. Although long tandem duplication is a very rare event in mtDNA duplication, the duplication can happen between the origin for H-strand replication (O H ) in the CR and the origin for L-strand replication (O L ) in the WANCY tRNA cluster, which is a distance of about twothirds of the genomic length. The mechanism for duplication between the CR and the WANO L CY tRNA genes in N. taihangnica could be caused by the O H and O L structures and be explained by the TDRL model [20], which is similar to the research of Shi et al. [24]. The hypothesized intermediate steps are as follows. Firstly, the above-mentioned O H and O L structures initiated DNA synthesis twice during mitochondrial replication, causing tandem duplication of the genes located between the CR and the WANCY region in the ancestral mitogenome (Fig. 5). Secondly, one of each of the duplicated gene pairs was randomly deleted completely or partially and then lost its function or became a pseudogene (Fig. 5).
In N. yunnanensis the genes of the mitogenome are extensively rearranged with clustering of eight genes on the L-strand in the same polarity and three control regions in an unexpected gene order. These special features of eight genes in the same polarity on the L-strand and two noncoding regions were reported in Crossorhombus azureus which proposed a new mechanism for gene rearrangement [23]. We can use this gene rearrangement mechanism to explain the polarity of gene rearrangement in N. yunnanensis. The hypothesized intermediate steps are as follows. Firstly, the inferred "dimer-mitogenome" intermediate of the N. yunnanensis mtDNA (Fig. 6) could be formed by two entire mitogenomes if the two mt genomes were linked by the head-totail method. Secondly, some duplicated genes were nonrandomly deleted completely except that all ten genes on the L-strand of one mt monomer were retained; some duplicated genes were also non-randomly deleted completely or partially from the other mt monomer (Fig. 6). Thirdly, the region of CR-trnP-trnQ-trnA-trnN-trnC-trnY-trnS-trnS-ND6-trnE was duplicated. Fourthly, some duplicated genes were randomly deleted completely. So the Dimer-Mitogenome and Non-Random Loss model (DMNR) [23] and the TDRL model [20] may be more appropriate to explain the gene arrangements in N. yunnanensis. But we have no suitable model to explain the four non-coding region (256 bp between trnQ and trnY, 722 bp between trnY and ND6, 1317 bp between trnP and trnA, 394 bp between trnC and trnS).

Phylogenetic analyses of Dicroglossidae
The evolutionary relationships of dicroglossid taxa indicated by the phylogenetic trees were mostly similar to previously reported molecular phylogeny [5]. Roelants et al. [78] suggested that Occidozygini is a sister clade to ((Dicroglossini + Paini) + Limnonectini), whereas van der Meijden et al. [79] found that Occidozyga (Occidozygini) is located within Dicroglossinae. Dubois [37] returned Occidozygini to Dicroglossinae as a tribe based on the strength of evidence produced by van der Meijden et al. [79]. In the present study, Occidozygini was found to be a sister clade to (Dicroglossini + (Paini + Limnonectini)), and Occidozygini (Occidozyginae) was observed to be a basal clade to Dicroglossinae.
In phylogenetic trees, the clade of (N. quadranus + ((N. taihangnica + N. taihangnica (KJ569109)) was clustered into the Nanorana. Although N. taihangnica and N. quadranus belong to the genus Feirana according to Fei et al. [42], we draw the conclusion that genus  The dimeric molecule with two monomers linked head-to-tail and subsequent first deletions of partial genes or complete genes resulting in the derived gene order. c State in N. yunnanensis after the first duplication and deletion. d Tandem duplication in the area from the control region (CR) to trnE (Glu) and subsequent second deletions of partial genes or complete genes resulting in the derived gene order. e State in N. yunnanensis after the second duplication and deletion. f Tandem duplication in the region of CR1-trnQ (Gln), CR2-trnP (Pro) and trnS (Ser)-CR3, respectively, and subsequent third deletions of partial genes or complete genes resulting in the derived gene order. g State in N. yunnanensis after the third duplication and deletion. * means that the tandem replicated genes were partially deleted and # means that the tandem replicated genes were completely deleted  Fig. 5 Proposed mechanism of gene rearrangements in N. taihangnica under a model of tandem duplication of gene regions and subsequent gene deletions. a Typical Dicroglossidae gene order. b Tandem duplication in the area from the control region (CR) to trnY (Tyr) and subsequent deletions of partial genes or complete genes resulting in the derived gene order. c State in N. taihangnica. * means that the tandem replicated genes were partially deleted and # means that the tandem replicated genes were completely deleted Feirana is not valid according to the phylogenetic relationship of Nanorana and Feirana, which was also supported by Frost et al. [25,43] and Che et al. [35,44] Invalidation of Q. robertingeri as a species The validity of Quasipaa robertingeri is also heatedly debated. Che et al. [35] found that Quasipaa robertingeri nested deeply within Q. boulengeri and suggested that Q. robertingeri should be synonymous with Q. boulengeri, which is supported by Frost et al. [25]. However, Fei et al. [36,42] insisted on the validity of Q. robertingeri as a species. The data of Pyron and Wiens [29] supported the proposal that Q. robertingeri was a sister clade to Q. shini, not to Q. boulengeri. To compare the genetic divergence we analyzed the complete mt genomes and 16S RNA gene of Q. boulengeri and Q. robertingeri in Mega 5.0 with the parameter p-distance model. The average genetic distance between Q. boulengeri and Q. robertingeri using mt genomes and 16S RNA was determined to be 4.3% and 1.1%, respectively, which is lower than the lowest interspecies mt genomes between Q. spinosa and Q. exilispinosa (6.8%) and 16S RNA diversity as a species threshold (3%) [80], respectively. Although Q. boulengeri is as a sister clade to Q. robertingeri in phylogenetic relationship, the genetic distance between Q. boulengeri and Q. robertingeri is lower than the genetic distance between interspecies of Quasipaa. The different gene arrangement of Q. boulengeri and Q. robertingeri cannot be used as a species delimitation method because the gene rearrangement can also happened within intraspecies. So we deduce that Q. robertingeri may not be a valid species.

Conclusion
The characteristics of mt genomes and gene arrangements provide novel insights into the phylogenetic relationships among several major lineages of Dicroglossidae. The phylogenetic relationship of ((Ranidae + Dicroglossidae) + (Mantellidae + Rhacophoridae)) is supported in BI analyses. Feirana is not a valid genus according to the phylogenetic relationship with Nanorana. Quasipaa robertingeri may be an invalid species according to genetic divergence. The gene arrangements of N. taihangnica and N. yunnanensis differed from those of other published dicroglossid mt genomes. The mt genomes are promising markers for discussing the reasons for intraspecies gene rearrangements, and the current results broadens our knowledge of the evolution of anuran mt genomes.