The first complete mitochondrial genome sequence of Nanorana parkeri and Nanorana ventripunctata (Amphibia: Anura: Dicroglossidae), with related phylogenetic analyses

Abstract Members of the Nanorana genus (family Dicroglossidae) are often referred to as excellent model species with which to study amphibian adaptations to extreme environments and also as excellent keystone taxa for providing insights into the evolution of the Dicroglossidae. However, a complete mitochondrial genome is currently only available for Nanorana pleskei. Thus, we analyzed the complete mitochondrial genomes of Nanorana parkeri and Nanorana ventripunctata to investigate their evolutionary relationships within Nanorana and their phylogenetic position in the family Dicroglossidae. Our results showed that the genomes of N. parkeri (17,837 bp) and N. ventripunctata (18,373 bp) encode 13 protein‐coding genes (PCGs), two ribosomal RNA genes, 23 transfer RNA (tRNA) genes, and a noncoding control region. Overall sequences and genome structure of the two species showed high degree of similarity with N. pleskei, although the motif structures and repeat sequences of the putative control region showed clear differences among these three Nanorana species. In addition, a tandem repeat of the tRNA‐Met gene was found located between the tRNA‐Gln and ND2 genes. On both the 5′ and 3′‐sides, the control region possessed distinct repeat regions; however, the CSB‐2 motif was not found in N. pleskei. Based on the nucleotide sequences of 13 PCGs, our phylogenetic analyses, using Bayesian inference and maximum‐likelihood methods, illustrate the taxonomic status of Nanorana with robust support showing that N. ventripunctata and N. pleskei are more closely related than they are to N. parkeri. In conclusion, our analyses provide a more robust and reliable perspective on the evolutionary history of Dicroglossidae than earlier analyses, which used only a single species (N. pleskei).

life history. However, due to their high-elevation habitats, Nanorana species experience extremely harsh abiotic factors, including hypoxia, high UV radiation, and dramatic temperature changes on a daily basis. Consequently, Nanorana is an excellent model species for studying the adaptations of frogs to extreme environmental conditions (Sun et al., 2015). Does the unique high-elevation environment of Nanorana have a greater impact on species differentiation and gene sequences characteristics? Our study aimed to clarify the mitochondrial genome sequence characteristics and phylogenetic relationship and the taxonomic status of the three species in the genus

Nanorana.
The phylogenetic relationships of Nanorana have been studied previously (Che et al., 2009;Chen, Wang, Liu, Xie, & Jiang, 2011;Lu, 1995;Zhou et al., 2014); however, debates on the taxonomic status of the three species that are the focus of this study are still ongoing.
The taxonomy of Nanorana species is not yet fully settled because of numerous changes during the last decade. Previous phylogenetic analyses support N. pleskei and Quasipaa spinosa as having a close relationship (Chen et al., 2011), as well as N. ventripunctata and N. parkeri (N. pleskei + (N. parkeri + N. ventripunctata)) (Lu, 1995). In other literature, however, N. pleskei and N. ventripunctata are reported to have a closer relationship (N. parkeri + (N. pleskei + N. ventripunctata)) (Che et al., 2009;Zhou et al., 2014), while Pyron and Wiens (2011) thought that N. pleskei and N. parkeri had a closer evolutionary relationship (N. ventripunctata + (N. pleskei + N. parkeri)). Thus, complete sequencing of the mtDNA in Nanorana can help clarify the phylogenetic relationships and genetic diversity within the genus. With that information, we can then better understand the phylogenetic status and intraspecific relationships among the three species within this group (Che et al., 2009;Chen et al., 2011;Jiang & Zhou, 2001Roelants, Jiang, & Bossuyt, 2004).
Mitochondrial genes such as the COX I, Cytochrome b (Cytb), D-loop, tRNA, and NADH have been used for previous phylogenetic and phylogeographic studies on the genetic divergence of Nanorana Liu et al., 2015;Wang et al., 2013;Zhang et al., 2010;Zhou et al., 2014). Here, we use complete mitochondrial genomes to analyze the phylogenetic relationships of the three Nanorana species (N. parkeri, N. ventripunctata, and N. pleskei) and other related species. Moreover, in order to reconstruct a robust evolutionary relationship among the three species, we need additional mitochondrial genomic information from Nanorana species. Therefore, we sequenced the complete mitochondrial genome of N. parkeri and N. ventripunctata and summarized the structural variations of 40 mitochondrial genome sequences in the Family Dicroglossidae. We reconstructed the phylogenetic relationships of Dicroglossidae using the concatenated sequences of 13 proteincoding genes from Dicroglossidae mitochondrial genomes, based on which the evolutionary characteristics of the mitochondrial genomes in Dicroglossidae were evaluated. Furthermore, we analyzed the mitochondrial genomic sequence and phylogenetic relationships within N. pleskei, N. ventripunctata and N. parkeri to assess the evolutionary status of the three species within the Nanorana genus.
Additionally, the complete mitochondrial genomes of two Nanorana species (N. ventripunctata and N. parkeri) were analyzed to find novel data with which to investigate the placement of the three Nanorana species in the phylogenetic tree of Dicroglossidae and to provide molecular data for further study on the taxonomic status and adaptive evolutionary mechanisms of these high-altitude species.

| Sampling and DNA extraction
The Xizang Plateau frog (N. parkeri, Figure 1) was sampled from Dangxiong County (4,300 m asl), the Tibet Autonomous Region, China, in September 2015. The Yunnan slow frog (N. ventripunctata) was sampled from Xianggelila County (4,200 m asl), Yunnan province, China, in July 2016. All collections were initially preserved in 95% ethanol and stored at −70°C until DNA extraction was performed. According to the protocol adopted by Zhang, Chao, Lai, Li, and Zhao (2000) and Xia, Liu, and Lu (2002), total mtDNA of two Nanorana species was extracted from skin tissue for the following PCR amplification.

| Mitochondrial DNA amplification and sequencing
The entire mitochondrial genome was amplified in twelve overlapping segments by PCR with LA-Taq DNA Polymerase (TaKaRa, China), using 20 ng of total genomic DNA from the sample as a template. Complete mtDNA was amplified as concatenated sequences by adopting selectively amplified mtDNA templates and 10 primer pairs, as published by Kurabayashi and Sumida (2009). Partial PCR primers were also designed based on the alignments of the relatively conserved regions of congeneric N. pleskei (NC_016119) and N. taihangnica (NC_024272). The PCR amplification was performed as follows: 2.5 min at 94°C, followed by 30 cycles of 0.5 min at 94°C, 0.5 min at 50-59°C, 3-5 min at 60°C, and a 9 min final extension at 72°C. PCR reactants were loaded on 0.8%-1.0% agarose gels, stained with ethidium bromide and photographed under ultraviolet light. PCR products were purified with Gel Extract Purification Kits (V-gene) and automated sequencing using an ABI 3730 sequencer, either directly or following subcloning into the pMD18-T vector (TaKaRa, China). To ensure maximum accuracy, each amplification product was sequenced twice independently, followed by a third PCR amplification.

| Sequence assembly and analysis
Sequences were assembled manually and aligned, and each gene was then translated into an amino acid sequence using MEGA 6.0 (Tamura, Stecher, Peterson, Filipski, & Kumar, 2013). The amino acid sequence alignments of each of the protein-coding genes (PCGs) were generated using the computer program Clustal X 1.83 (Thompson, Gibson, Plewniak, Jeanmougin, & Higgins, 1997). Based on sequence similarity results from BLAST searches, ribosomal RNA (rRNA) genes were recognized, and tRNA genes were identified using tRNAscan-SE 1.21 (Schattner, Brooks, & Lowe, 2005). Base composition and codon usage were analyzed in MEGA 6.0 (Tamura et al., 2013). The mitochondrial genome sequences have been submitted to NCBI GenBank with the accession number NC_026789 (N. parkeri) and KY594708 (N. ventripunctata). Features of the base composition of nucleotide sequences were detected using the ATskew and GC-skew in the mitochondrial genome. We then calculated the AT-skew and GC-skew using the following formulae from Perna and Kocher (1995): AT-skew = (A − T)/(A + T) and GC-skew = (G − C)/ (G + C).

| Phylogenetic analysis
Combined with 38 other Dicroglossidae mitochondrial genomes from NCBI GenBank (Supporting information Table S1), the mitochondrial genomes of N. parkeri and N. ventripunctata were analyzed using the phylogenetic tree method, with the concatenated sequences of the 13 protein-coding genes and the two species Babina subaspera (NC_022871) and Hylarana guentheri (NC_024748) as outgroups. First, we aligned the 13 mitochondrial protein-coding gene sequences in Clustal X 1.83 (Schattner et al., 2005) with the default settings, and then we concatenated individual genes excluding the stop codon. We selected the optimal nucleotide substitution model in jModeltest v0.1.1 (Posada, 2008) and used the Akaike Information Criterion (AIC: Posada & Buckley, 2004).
Maximum likelihood (ML) and Bayesian inference (BI) were used for phylogenetic analyses in MrBayes 3.2.2 (Ronquist et al., 2012), and BI of nucleotide acid datasets was performed using the GTR + I + G model (Lanave, Preparata, Saccone, & Serio, 1984). A ML tree was constructed using RAxML, and the robustness of the phylogenetic results was tested through bootstrap analysis with 1,000 replicates (Stamatakis, 2014).
The gene components were very loosely juxtaposed with 134/42 (N. parkeri) and 63/39 (N. ventripunctata) of gap/overlapping nucleotides, compared to that of N. pleskei (71/49; Table 1) (Simon et al., 1994). Although the overall A + T contents of 57.87% in N. parkeri and 59.1% in N. ventripunctata were relatively higher than that of N. pleskei (57.5%), those values are within the range (52.8%-62.74%) of Dicroglossidae (Supporting information Table S1). The nucleotide skew was highly similar among these mitochondrial genomes including that of N. pleskei, with only some exceptions found on COX2, ATP6, ATP8 and the putative control region (Table 2).

| Protein-coding genes (PCGs) and codon usage patterns
The inferred start/stop codons for protein-coding genes of N. parkeri, N. ventripunctata, and N. pleskei are listed in Table 1. In three mitochondrial genomes, the protein-coding genes were initiated by ATG, with the exceptions of COX1, ND1, ND2, and ND3 (Table 1). The open reading frame of ND1 and ND3 started with GTG, while that of COX1 and ND2 started with ATA and ATT, respectively. The canonical stop codon (TAA or TAG) can be found in four protein-coding genes (ATP8, Cytb, ND4L, and ND5; Table 1), while COX1 and ND6 use AGG and AGA as the termination codon, respectively. The remaining seven (ATP6, COX2-3, and ND1-4) had incomplete T-stop codons (Table 1), completed (TAA) by polyadenylation after transcription (Boore, 2001).
The relative synonymous codon usage (RSCU) values of the three species of Nanorana mitogenomes are shown in Table 3, Supporting information Tables S2-S4. The results demonstrate that synonymous codon usage has a distinct bias toward A or T for 13 PCGs. The codons AUU (5.03%-5.62%), UUU (3.92%-4.37%), GCC (3.84%-4.03%), and CUU (3.63%-3.79%) were the four most frequently used codons in the mitogenomes of our three species of Nanorana, accounting for 16.42%-17.81%. In addition to GCC codon, these codons were mainly composed of A or U nucleotides, indicating the highly biased usage of A and T nucleotides in the three species of Nanorana PCGs. Meanwhile, the most frequently represented amino acids in the three species of Nanorana mitochondrial proteins were Leu (16.27%-16.38%), Ala (8.24%-8.32%), Ile (7.95%-8.11%), and Phe (6.62%-6.78%), accounting for 39.08%-39.59%. The least frequently represented amino acid was Cys (0.74%-0.77%). Codon usage of PCGs showed a major bias of A + T content, which played a major role in the A + T bias of the entire mitogenome. Similar patterns with a strong T-or A-bias in the wobble position have been found among other Nanorana species also. The RSCU analysis showed that codons with A or T (U) at the third position are mostly overused compared with other synonymous codons. Therefore, the codon usage can reveal nucleotide bias too. These data imply a high A + T content in the three Nanorana species. The bias toward the use of Ts over As, to the 13 PCGs, is more obvious in these three Nanorana mitogenomes with −0.080 to −0.100 AT skewness. Moreover, negative AT-skew and GC-skew were found in the third position, whereas both the first and second positions showed positive AT-skew and negative GC-skew in N. parkeri and N. ventripunctata. In contrast, the first, second and third positions showed negative AT-skew and GC-skew in N. pleskei (Supporting information Table S5). genes. Extra tRNA-Met was also found in Quasipaa boulengeri, Fejervarya cancrivora, Hoplobatrachus rugulosus, Euphlyctis hexadactylus, Limnonectes bannaensis, and Occidozyga martensii (Alam et al., 2010;Chen et al., 2011;Li et al., 2014aLi et al., , 2014bRen et al., 2009;Shan, Xia, Zheng, Zou, & Zeng, 2014;Zhang et al., 2009).

| Transfer and ribosomal RNA genes
But this phenomenon is different to that seen in Amolops tormotus and other typical vertebrates (Su, Wu, Yan, Cao, & Hu, 2007).
Two tRNA-Met genes in each lineage may come from different origins (Kurabayashi et al., 2006), and the tandem duplication of the tRNA-Met gene can be seen as a synapomorphic feature of Dicroglossidae. A tandem duplication of the mitochondrial tRNA-Pro and tRNA-Thr genes in Bipes biporus has been reported from previous research (Macey, Schulte, Larson, & Papenfuss, 1998).

| Noncoding regions
Putative control region, of 2,259 bp (N. parkeri) and 2,857 bp (N. ventripunctata) were found in Cytb and tRNA-Leu, which is longer than that of N. pleskei (2,143 bp) ( Table 1). The size of control region variation demonstrated different lengths of the total mitogenomes for the three species. The A + T contents (65.96% in N. parkeri and 69.86% in N. ventripunctata) in control region were higher than in other regions (Table 1). Additionally, the A + T contents rated different lengths of the total mitogenomes for the three species. The A + T content in this region is higher than that in the coding regions (Boore, 1999;Simon et al., 1994).    (Table 4).
Three conserved sequence blocks (CSBs) may be related to in the initiation of the mtDNA synthesis and they (CSB-1, CSB-2, CSB-3) can be identified between the tandem repeat units at the 5′ and 3′-sides (Table 4; Figure 4). CSB-1, CSB-2 and CSB-3 of N. ventripunctata and N. parkeri showed high similarity to the consensus in other amphibians, while the variation in N. pleskei is slightly larger ( Figure 5); moreover, CSB-1 is not reduced to a truncated penta motif (5′-GACAT-3′) as it is in the caecilians (San Mauro et al., 2004;Zardoya & Meyer, 2000). However, a truncated CSB-1 had been recorded in Xenopus laevis (Anura) (Roe, Ma, Wilson, & Wong, 1985). The CSB-2 motif was not found in N. pleskei ( Figure 5). In addition, the multiple motifs of mtDNA control regions (CR) may be associated with the transcription and replication of the mitochondrial genome (Taanman, 1999). The function of these conserved sequence blocks is unclear. Further study on the mechanistic basis of mtDNA replication is warranted for Nanorana species.

| Phylogenetic relationships
The concatenated PCG data of the mitogenome sequences in our study contained 11,292 nucleotide positions, including 4,314 conserved sites, 6,978 variable sites and 6,505 potentially parsimony-informative sites. Phylogenetic trees were reconstructed using BI and ML analyses, based on the nucleotide dataset. The use of PCG sequences of the mitogenomes has become an informative strategy for inferring phylogenetic relationships (Boore, Macey, & Medina, 2005). Using the 13 PCG sequences to concatenate may achieve a more complete analysis. BI and ML methods consistently support similar tree topologies by strong node-supporting values.
So far, combined with the 38 mitochondrial genome sequences in GenBank database, our phylogenetic analyses revealed that the subfamily Dicroglossinae's monophyly was well supported (Li et al., 2014a(Li et al., , 2014bRoelants et al., 2007;Yuan et al., 2016). The subfamily Dicroglossinae is the sister clade to the Occidozyginae ( Figure 6).
The Dicroglossinae species was divided into two clades with one clade (Clade 1) containing Nanorana, Quasipaa, and Limnonectes, and the other (Clade 2) including Fejervarya, Euphlyctis, and Hoplobatrachus ( Figure 6), as supported by previous studies (Lv, Bi, & Fu, 2014;Yuan et al., 2016;Zhang, Xia, & Zeng, 2016). Quasipaa and Nanorana be-  Table 4 the Nanorana species clustered in a single monophyly. Our molecular phylogeny indicates N. ventripunctata and N. pleskei are more closely related compared with N. parkeri, and strongly supports that N. parkeri is basal to N. pleskei and N. ventripunctata based on 13 PCG genes of the mitogenome (BP = 100%, PP = 1.00) (Figure 6), in agreement with the relationships inferred by the research report of Che et al. (2009Che et al. ( , 2010 contain species-specific evolutionary affinities, which can be efficiently recovered by improving taxon sampling (Rubinstein et al., 2013).

| CON CLUS IONS
In summary, the complete mitochondrial genomes of two Nanorana species were determined in this study. Our mitogenome analyses, including gene content, gene order, strand asymmetry, base composition, rRNA and tRNA secondary structure and phylogenetic analysis, indicate several significant features: a tandem repeat of the tRNA-Met gene was detected in three Nanorana species, located between the tRNA-Gln and ND2 genes. The control region contains distinct repeat regions at both 5′ and 3′-sides, and the CSB-2 motif was not found in the N. pleskei. Based on nucleotide sequences of 13 PCGs, and using BI and maximum-likelihood F I G U R E 6 Results of phylogenetic analyses using BI and ML analysis indicated evolutionary relationships among 38 individuals based on 13 PCGs sequences. Babina subaspera (NC_022871) and Hylarana guentheri (NC_024748) were used as outgroups. Tree topologies produced by BI and ML analyses were equivalent. Bayesian posterior probability (PP) and bootstrap support (BP) values for ML analyses are shown in order on the nodes. The asterisks indicate new sequences generated in this study F I G U R E 5 Structures and alignments of identified putative terminationassociated sequences (TAS) and, conserved sequence blocks (CSB 1-3). Alignment gaps and nucleotides identical to the first line are indicated by dashes (-) and a dot (•), respectively. Variable nucleotides are marked with corresponding nucleotides analyses, the phylogenetic data illustrate the taxonomic status of Nanorana and provides robust support that N. ventripunctata and N. pleskei are more closely related than N. parkeri. Our study provides useful additional data for further phylogenetic analysis of the Nanorana genus. Expanding our knowledge of the phylogenetic relationships within the Nanorana genus will ultimately aid in future research to protect and maintain biodiversity within many other anuran species. However, the proposed evolutionary relationships among these three species based on the findings that emerged in the study should be accepted with caution due to limited taxon sampling. Many aspects of the phylogeny of the genus Nanorana remain to be resolved and further analysis based on more molecular information (including nuclear gene data) and extensive taxon sampling is necessary to elucidate the phylogenetic relationships among genus Nanorana or Dicroglossidae.

ACK N OWLED G M ENTS
We thank Xiaoyan Ma for providing the photograph of Nanorana parkeri and two anonymous reviewers for helpful comments that

CO N FLI C T O F I NTE R E S T
The authors have declared that no competing interest exists.

AUTH O R CO NTR I B UTI O N
Jiang L., Ruan Q., and Chen W. designed the manuscript, You Z. and Yu P. analyzed the data, and Jiang L. and Chen W. wrote the manuscript.