Comparative Mitogenomics of True Frogs (Ranidae, Anura), and Its Implications for the Phylogeny and Evolutionary History of Rana

Simple Summary The true frogs of the genus Rana are a complex and diverse group. Many new species have been discovered with the help of molecular markers and morphological traits. However, the evolutionary history in Rana were not well understood. In this study, we sequenced and annotated the complete mitochondrial genome of R. longicrus and R. zhenhaiensis. In 13 protein codon genes, the COI was the most conserved, and ATP8 had a fast rate of evolution. The Ka/Ks ratio analysis among Rana indicated the protein-coding genes were suffering purify selection. There were three kinds of gene arrangement patterns found. This study provides mitochondrial genetic information, improving our understanding of mitogenomic structure and evolution, and recognizes the phylogenetic relationship and taxonomy among Rana. Abstract The true frogs of the genus Rana are a complex and diverse group, containing approximately 60 species with wide distribution across Eurasia and the Americas. Recently, many new species have been discovered with the help of molecular markers and morphological traits. However, the evolutionary history in Rana was not well understood and might be limited by the absence of mitogenome information. In this study, we sequenced and annotated the complete mitochondrial genome of R. longicrus and R. zhenhaiensis, containing 22 tRNAs, 13 protein-coding genes, two ribosomal RNAs, and a non-coding region, with 17,502 bp and 18,006 bp in length, respectively. In 13 protein codon genes, the COI was the most conserved, and ATP8 had a fast rate of evolution. The Ka/Ks ratio analysis among Rana indicated the protein-coding genes were suffering purify selection. There were three kinds of gene arrangement patterns found. The mitochondrial gene arrangement was not related to species diversification, and several independent shifts happened in evolutionary history. Climate fluctuation and environmental change may have played an essential role in species diversification in Rana. This study provides mitochondrial genetic information, improving our understanding of mitogenomic structure and evolution, and recognizes the phylogenetic relationship and taxonomy among Rana.


Introduction
The Ranidae (Anura) is comprised of 24 genera and approximately 434 species worldwide [1]. The true frogs of the genus Rana are a complex and diverse group that are widely In this study, we sequenced and annotated the complete mitogenome of two Rana species (R. longicrus and R. zhenhaiensis). The genome composition and characteristics were analyzed, the relative synonymous codon usage (RSCU) and AT-skew values of the PCGs were calculated. Furthermore, we analyzed phylogenetic relationships, gene arrangement, and selective pressures among Rana species. In the future, sequencing more mitogenomes from various taxonomic levels will provide useful information for better understanding the evolutionary history and will contribute to further taxonomic research within Ranidae.

Sample Collection and DNA Extraction
The samples of R. longicrus was collected from Chiu-ling Mountains, Jiangxi Province, China (29 • 56.64 E). After collection, the muscles were initially preserved in 75% ethanol for two days, then preserved in 95% ethanol, and finally transferred to −20 • C in the laboratory for long-term storage in Nanjing Normal University. Total genomic DNA was extracted using a DNAeasy tissue kit (Qiagen, Hilden, Germany) following the manufacturer's instructions.

Library Preparation and Sequencing
The extracted DNA was then sent to Novogene (Beijing, China) for sequencing on the Illumina sequencing platform. The concentration of DNA was checked with a Nanodrop 1000 Spectrophotometer (Thermo Scientific, Waltham, MA, USA). The extracted DNA was sheared to 400-600 bp using an ultrasonic technique. The sequencing library was produced using the Illumina Truseq DNA Sample Preparation Kit (Illumina, San Diego, CA, USA) according to the manufacturer's instructions. The prepared library was loaded on the Illumina Novaseq 6000 platform for PE 2 × 150 bp sequencing at Novogene (Beijing, China).
A total of 6.0 Gb raw reads were generated by next-generation sequencing on the Illumina platform. Before assembly, Illumina raw data were filtered into clean reads, and then undesirable reads were removed by fastp v. 0.21 [23], with the following parameters "-q 15 -u 40 -5 -x -w 40 -f 10 -F 10" [23]. This filtering step was performed to remove the reads with adaptors, the reads showing a quality score below 20 (Q < 20), containing a percentage of uncalled based ("N" characters) equal or greater than 10% and the duplicated sequences. After removing low quality sequences, the 5.67 Gb clean Raw Data of R. zhenhaiensis and 5.05 Gb of R. longicrus were harvested. The clean reads de novo assemblies were conducted in Geneious 10.1.2 using the mitogenome of R. dybowskii (GenBank no. KF898355) as a reference map [24], and aligned contigs (≥80% similarity and an average 170 X coverage) were ordered according to the reference genome. The assembled contig identified as mitogenome was manually examined for repeats at the beginning and end of the sequence to infer circularity.

Sequence Assembly and Annotation
The mitogenome was first annotated with MITOS [25]. The PCGs and rRNA genes were determined by NCBI open reading frame finder implemented at the NCBI website with the vertebrate mitochondrial genetic code, and then manually corrected by comparison with the available sequences of closely related Rana species downloaded from GenBank using ClustalW in MEGA X [26]. The, tRNAscan-SE [27], and ARWEN [28] were used to confirm tRNA and rRNA annotation results. All tRNA genes were identified by their cloverleaf secondary structure using tRNA-scan SE or determined by comparison with the homologous sequences [27]. The graphical map of the mitogenome map was conducted using CGView Server online software [29].
The software of MEGA X was used to calculate the number of variable sites, the parsimony informative sites, the singleton, and the average uncorrected pairwise distances for 13 protein-coding genes of Rana [26]. Codon usage was estimated using DnaSP 5.1 [30]. Nucleotide composition and the relative synonymous codon usage (RSCU) were calcu-lated with MEGA X [26]. The rates of non-synonymous substitutions (Ka, π modified), synonymous substitutions (Ks, π modified), the effective number of codons (ENC) and the codon bias index (CBI) for each protein-coding gene were determined with DnaSP 6.0 [31]. The value of nucleotide composition skewness was measured using the following formulas: GC-skew = (G − C)/(G + C) and AT-skew = (A − T)/(A + T) [32,33]. The tandem repeats were searched in the CR using the Tandem Repeats Finder program (https://tandem.bu.edu/trf/trf.html, accessed on 18 November 2021) [34].

Molecular Phylogenetic Analysis
We constructed the phylogenetic topology of 18 Rana species using 13 PCGs, 12S and 16S with two species (Babina subaspera GenBank no. NC_02287 and Pelophylax nigromaculatus KT878718) as outgroups. Phylogenetic analysis was performed using Bayesian Inference (BI). To determine the optimal partitioning of the data, the best-fit partitioning scheme and the most appropriate nucleotide evolution model for each partition were implemented in PartitionFinder 2, with greedy algorithm and Akaike Information Criterion (AICc) criterion [35]. BI method was performed using MrBayes 3.1.2 [36]. Four Markov Chains Monte Carlo (MCMC) chains were run for 1.0 × 10 6 generations. Two independent runs were performed to confirm consistent approximation of the posterior parameter distributions. Stationarity was reached when the average standard deviation of split frequencies was below 0.01. The convergence of MCMC runs and effective sample sizes (ESS > 200) were checked by plotting the log-likelihood scores against the generation times using the program TRACER 1.6 (http://beast.bio.ed.ac.uk/Tracer accessed on 18 November 2021). The first 25% of sampled trees and estimated parameters were discarded as burn-in. The remaining trees were used to calculate consensus tree and Bayesian posterior probabilities.

Divergence Times Estimates
Divergence times were estimated using BEAST 1.8.4 [37]. The following parameters were used: GTR + I + G substitution model, relaxed log-normal molecular clock model and a Yule process for tree prior. Two calibration points were used for time calibration. The first calibration point, representing the split between Eurasian and American species, was set as 31.2 ± 8.1 Ma [2,38]. The second calibration point was a constraint of 15.0 Ma with 2.5% and 97.5% quantiles of 15.2 and 23.6 Ma for the most recent common ancestor of the R. catesbeiana group [2,13,39]. We performed two independent runs with a MCMC chains of 2 × 10 8 generations, and trees were sampled every 1000 generations. The first 25% of the generations were discarded as burn-in. Convergence of the chains was determined with TRACER ver. 1.7.1 [40], with target ESS values more than 200 for all parameters. The tree information was annotated and visualized by FigTree v. 1.4.4 (http://tree.bio.ed.ac.uk/ software/figtree/, accessed on 24 October 2021).

Mitogenomic Rearrangements Analysis
Multiple sequence alignment of the Rana mitogenome sequences was performed to verify the gene order of the mitochondrial genome. We reconstructed the ancestral states of gene rearrangements. The rearrangements patterns were categorized into three types (result 3.4). Ancestral discrete characteristics were reconstructed in a time-calibrated MCMC tree using the "ace" function in the "ape" package [41]. The parameters of the models were estimated using the ML method [42]. The final states were plotted in the time-calibrated MCMC tree using the "phytools" package [43].

General Characteristics of the Mitogenome
The complete mitogenome of R. longicrus and R. zhenhaiensis was 17,502 bp and 18,006 bp in length, respectively. They carried the typical composition, including 37 genes (13 PCGs, two rRNAs, and 22 tRNAs) and the control region ( Figure 1). Length differences were primarily the result of variation in control region. The overall base composition of the

Protein-Coding Genes and the Codon Usage
Relative synonymous codon usage (RSCU) values for the 13 PCGs were summarized in Table 2. The 13 most frequently used amino acid was Ile AUU , which accounts for 4.56% of the usage, and the least was Arg CGG (0.17%). To further study the codon usage bias of Rana, the correlation between ENC (effective codon number), CBI (codon bias index), G + C content of all codons (G + Cc), and G + C content of the third codon position (G + C3s) were analyzed. We found a significant negative correlation between CBI and ENC (R 2 = 0.71; p < 0.01). However, there was no correlation between the other pairs ( Figure 2). The total length of the protein-coding genes in each species was 11,211 bp after removing termination codons and indels. The length of 13 PCGs ranged from 159 bp (ATP8) to 1740 bp (ND5). The overall AT-skew and GC-skew of the 13 PCGs were negative except that the AT-skew of COII and ATP8 and the GC-skew of ND6 were positive ( Table 3). The use of start codons for 13 PCGs were quite common. Four start codons (ATG, GTG, ATC and ATA) were detected in 13 PCGs. The most common start codon was ATG, which accounts for 70.94% of the start codons, followed by GTG (16.67%) ( Figure 3). However, the start codon is not determined in ND1 in R. longicrus and R. zhenhaiensis. Four types of stop codons were complete stop codons (TAA, TAG, AGA, and AGG), and the other type was incomplete stop codons (T-) ( Figure 3). The incomplete stop codons were the most common stop codon, which was used by 10 genes (ND1, ND2, COII, ATP8, ATP6, COIII, ND3, ND4, ND5, Cyt b).      Comparison of PCGs provide a better understanding of the evolutionary pattern of molecular evolution. The variable positions in each gene ranged from 35.50% (COIII) to 54.09% (ATP8), and the parsimony informative sites ranged from 30.27% (COIII) to 47.41% (ND5). The average uncorrected pairwise distances showed the heterogeneity of the evolutionary rate in each gene. The COI (0.02), COII (0.05), and COIII (0.05) had slowly evolutionary rate, whereas the gene of ATP8 (0.21) had the fastest rate (Table 3).
To further understand the evolutionary patterns among the 13 protein-coding genes, Ka/Ks was calculated for each protein-coding gene, respectively. The Ka/Ks ratio values of 13 PCGs were <1, which suggested purifying selection of these functional genes. The gene of COI (0.02) had the slowest evolution rates, and was then followed by COIII (0.04) and COII (0.04). The gene of ATP8 (0.18) had the fastest evolutionary rates and was then followed by ND3 (0.17), ND4L (0.15), and ND6 (0.14) ( Table 3).

Transfer RNA, Ribosomal RNA Genes and Control Region
The mitogenome of R. longicrus and R. zhenhaiensis includes 22 tRNA genes, with the length ranging from 65 to 73 bp. The majority of tRNAs exhibit a secondary structure, usually a cloverleaf shape, with alternating double-helix stems and single-stranded nucleotides. However, tRNA Ser(AGY) absented the dihydrouridine (DHU) arm.
The control region is an important non-coding region in the mitogenome. It is located between Cyt b and tRNA Leu(CUN) , which spans 1,879 bp and 2,346 bp with A + T content of 61.6% and 59.3%, respectively. In R. zhenhaiensis, only one tandem repeats sequence of motif-1: 5 -TATGTTTAATAATCATTAATCTATCTGGATACTATCTC-3 (38 nucleotides with 5 tandem repeats) was found. In R. longicrus, two tandem repeats were found, the motif-1 had 38 nucleotides with 4 tandem repeats, and the motif-2: 5 -TATGTTTAATAATCATTAACCTATCTAAGTACTATACCTA-3 had 40 nucleotides with 2 tandem repeats).

Phylogenetic Relationships
We conducted a phylogenetic analysis of the mitogenomes of Rana species by using BI methods, estimating the phylogenetic tree using nucleotide sequences of 13 PCGs, 12S and 16S. All of the nodes exhibited high posterior probability values (Figure 4). In the phylogenetic tree, there were two major clades among Rana. The Clade A was consisted of three species, with R. sylvatica forms the sister group of two frogs (R. catesbeiana and R. okaloosae). The Clade B was consisted of the remaining species. The species of R. longicrus and R. zhenhaiensis formed a clade, which was a sister group to R. omeimontis and R. chaochiaoensis. Our analyses resolved two geographic lineages of Rana: Clade A (North America) was the sister group of Clade B (East Asia). The phylogeny of Rana reconstructed in our study was very similar to previous studies [2,44]. The main conflict involved the relationships among three species (R.chensinensis, R. kukunoris and R. huanrensis). Yuan et al. (2016) revealed the relationship of ((R.chensinensis, R. huanrensis), R. kukunoris) with relative low nodal support using sequences of six nuclear and three mitochondrial loci (total of 7,250 bp) [2]. However, this study supported the relationship of ((R.chensinensis, R. kukunoris), R. huanrensis) with high nodal support value, which was in agreement with the result of Yang et al. (2018) using two rRNAs and 13 PCGs [14]. Due to limited data on Europe and Central Asia groups in this study, we cannot provide more effective results to resolve detailed phylogenetic problems. The newly sequenced mitogenomes will provide useful information for better understanding the evolutionary history of genus Rana and will contribute to further taxonomic research within Ranidae.
(total of 7,250 bp) [2]. However, this study supported the relationship of ((R.chensinensis, R. kukunoris), R. huanrensis) with high nodal support value, which was in agreement with the result of Yang et al. (2018) using two rRNAs and 13 PCGs [14]. Due to limited data on Europe and Central Asia groups in this study, we cannot provide more effective results to resolve detailed phylogenetic problems. The newly sequenced mitogenomes will provide useful information for better understanding the evolutionary history of genus Rana and will contribute to further taxonomic research within Ranidae.

Divergence-Time Estimation
The high ESS (effective sample size) value (>200) was identified for all parameters in the BEAST analysis to estimate the divergence time of Rana species ( Figure 5). Dating analyses suggests that the most recent common ancestor (TMRCA) of Rana dates back to at 28.17 million years ago (Mya) [mean value; 95% of the highest posterior density (HPD), 23.07-34.69 Mya] (95% HPD, Mya) and dated to Oligocene ( Figure 5), roughly consistent with the previous study [2]. It may be affected by the collision of India with Eurasia between 35-20 Ma [38]. TMRCA of Clade A and B were estimated at 22.00 Mya (95% HPD, 13.71-30.43 Mya) and 23.16 Mya (95% HPD, 19.19-328.11 Mya), respectively. Since 16.09 Mya, the diversification rate within Clade B began to increase during the Miocene ( Figure 5), which indicated that the ancestral Rana fast dispersal across East Asia [38]. This study inferred the separation of R. dybowskii and R. uenoi dated at 8.13 Mya (95% HPD, 5.41-11.11 Mya) during the Middle and Late Miocene, which was consistent with the previous study [13], and was similar to the species diversification pattern in Hyla [45]. The Early Miocene diversification of R. zhenhaiensis and R. longicrus was formed by volcanism in Sikhote-Alin and Sakhalin in mid-Miocene-Pliocen [46]. We suggest that climate fluctuation and environmental change may have played an essential role in the species diversification in Rana.
inferred the separation of R. dybowskii and R. uenoi dated at 8.13 Mya (95% HPD, 5.41-11.11 Mya) during the Middle and Late Miocene, which was consistent with the previous study [13], and was similar to the species diversification pattern in Hyla [45]. The Early Miocene diversification of R. zhenhaiensis and R. longicrus was formed by volcanism in Sikhote-Alin and Sakhalin in mid-Miocene-Pliocen [46]. We suggest that climate fluctuation and environmental change may have played an essential role in the species diversification in Rana.

Mitogenomic Rearrangements
Based on the comparison of genome organization, we detected three kind of gene arrangements. Our results showed that Pattern 1 was the most common type in Rana mitogenomes ( Figure 6). Reconstructions of mitogenomic rearrangements pattern indicated that the ancestor of Rana consisted of Pattern 1. The R. kunyuensis and R. amurensis species have the same rearrangement, indicating that they have close relationships. The results showed that the mitochondrial gene arrangement was not related to species diversification, and several independent shifts happened in evolutionary history. mitogenomes ( Figure 6). Reconstructions of mitogenomic rearrangements pattern indicated that the ancestor of Rana consisted of Pattern 1. The R. kunyuensis and R. amurensis species have the same rearrangement, indicating that they have close relationships. The results showed that the mitochondrial gene arrangement was not related to species diversification, and several independent shifts happened in evolutionary history. All rearrangements occurred between ND4 and tRNA Thr region (Figure 7). Our results showed that Pattern 1 (the typical Neobatrachian-type arrangement) was the most common type, and Pattern 2 was shared by R. kunyuensis and R. amurensis. Compared with Pattern 1, one additional D-loop region was inserted into the upstream of tRNA Thr , and the ND5 was translocated from the typical tRNA Ser(AGY) downstream to the tRNA LeuCUN downstream. Pattern 3 was only discovered in R. pyrenaica, in which tRNA Pro -tRNA Phe was not located between tRNA Thr and 12S.
The gene arrangement of vertebrate mitogenome is usually conservative, and gene recombination is relatively rare or random [49]. However, in Neobatrachian, tRNA genes were translocated from typical positions, and these tRNA genes formed a cluster upstream of the 12S rRNA gene [50]. Currently, the duplication and random loss (TDRL) model can be used to explain most of the animal mitogenome reorganization [51]. In this model, a part of the entire genome was duplicated accidentally due to replication errors (either slipped strand mispairing or inaccurate termination). Then one of the duplicates All rearrangements occurred between ND4 and tRNA Thr region (Figure 7). Our results showed that Pattern 1 (the typical Neobatrachian-type arrangement) was the most common type, and Pattern 2 was shared by R. kunyuensis and R. amurensis. Compared with Pattern 1, one additional D-loop region was inserted into the upstream of tRNA Thr , and the ND5 was translocated from the typical tRNA Ser(AGY) downstream to the tRNA LeuCUN downstream. Pattern 3 was only discovered in R. pyrenaica, in which tRNA Pro -tRNA Phe was not located between tRNA Thr and 12S.
The gene arrangement of vertebrate mitogenome is usually conservative, and gene recombination is relatively rare or random [49]. However, in Neobatrachian, tRNA genes were translocated from typical positions, and these tRNA genes formed a cluster upstream of the 12S rRNA gene [50]. Currently, the duplication and random loss (TDRL) model can be used to explain most of the animal mitogenome reorganization [51]. In this model, a part of the entire genome was duplicated accidentally due to replication errors (either slipped strand mispairing or inaccurate termination). Then one of the duplicates (or CR) was converted into a pseudogene and subsequently excised from the genome through the accumulation of natural mutations [52]. In this study, R. longicrus and R. zhenhaiensis mitogenome possessed the same gene order as Neobatrachian. Nine genes (between tRNA Leu (CUN) and tRNA Phe ) were tandem duplicated, forming two identical gene clusters. Then, after the random loss of nine genes, the gene cluster produced by tandem duplication formed the currently known type of mitochondrial gene rearrangement Figure 7b. This explanation was confirmed by previous studies [53,54].
(or CR) was converted into a pseudogene and subsequently excised from the genome through the accumulation of natural mutations [52]. In this study, R. longicrus and R. zhenhaiensis mitogenome possessed the same gene order as Neobatrachian. Nine genes (between tRNA Leu (CUN) and tRNA Phe ) were tandem duplicated, forming two identical gene clusters. Then, after the random loss of nine genes, the gene cluster produced by tandem duplication formed the currently known type of mitochondrial gene rearrangement Figure 7b. This explanation was confirmed by previous studies [53,54]. Various types of mitogenome recombination have occurred in Ranidae species, which is also reflected in Rana genus [49]. In this study, we found a total of three kind of gene arrangement patterns in Rana species, which was consistent with the various types of mitogenome recombination in Ranidae species [55]. The mitogenome of R. amurensis and R. kunyuensis possessed two duplicate D-loop regions and this phenomenon was also found in other frogs [53]. The duplicated D-loop regions were similar to the original Dloop structure, which will result from homologous recombination between paralogous Dloop regions [44]. In general, the observed pattern is largely consistent with the TDRL model, the most general explanation for mitogenomic reorganization [56]. Subsequently, gene transfers via retrotransposition may be a pattern of gene rearrangement in animal mtDNAs [55]. We found mitogenome rearrangement in genus Rana, but it remains to be seen whether similar or new rearrangement patterns will occur in un-sequenced Rana.

Conclusions
In this study, we sequenced and annotated the mitogenome of R. longicrus and R. zhenhaiensis using the Illumina Novaseq 6000 platform for PE 2 × 150 bp sequencing. Circular mitogenomes of Rana displayed high size variation, with a mean length of 18,534 bp, ranging from 16,961 to 22,255 bp. The length of the mitogenome of R. longicrus was smaller than that of most Rana species. Differences in the length of the complete mitogenome are mainly caused by the control region. The mitogenome encoded a control region and a typical set of 37 genes containing 2 rRNA genes, 13 protein-coding genes, and 22 tRNA genes. The genome composition and characteristics were analyzed, RSCU and AT-skew values of the PCGs were calculated. In general, negative AT-skew and negative GC-skew are found in Rana mitogenomes, implying the specific bias toward T Various types of mitogenome recombination have occurred in Ranidae species, which is also reflected in Rana genus [49]. In this study, we found a total of three kind of gene arrangement patterns in Rana species, which was consistent with the various types of mitogenome recombination in Ranidae species [55]. The mitogenome of R. amurensis and R. kunyuensis possessed two duplicate D-loop regions and this phenomenon was also found in other frogs [53]. The duplicated D-loop regions were similar to the original D-loop structure, which will result from homologous recombination between paralogous D-loop regions [44]. In general, the observed pattern is largely consistent with the TDRL model, the most general explanation for mitogenomic reorganization [56]. Subsequently, gene transfers via retrotransposition may be a pattern of gene rearrangement in animal mtDNAs [55]. We found mitogenome rearrangement in genus Rana, but it remains to be seen whether similar or new rearrangement patterns will occur in un-sequenced Rana.

Conclusions
In this study, we sequenced and annotated the mitogenome of R. longicrus and R. zhenhaiensis using the Illumina Novaseq 6000 platform for PE 2 × 150 bp sequencing. Circular mitogenomes of Rana displayed high size variation, with a mean length of 18,534 bp, ranging from 16,961 to 22,255 bp. The length of the mitogenome of R. longicrus was smaller than that of most Rana species. Differences in the length of the complete mitogenome are mainly caused by the control region. The mitogenome encoded a control region and a typical set of 37 genes containing 2 rRNA genes, 13 protein-coding genes, and 22 tRNA genes. The genome composition and characteristics were analyzed, RSCU and AT-skew values of the PCGs were calculated. In general, negative AT-skew and negative GC-skew are found in Rana mitogenomes, implying the specific bias toward T and C in nucleotide composition. The average uncorrected pairwise distances showed that the COI may be the most conserved protein coding gene, and ATP8 was the least conserved. The most common start codon was ATG and the incomplete stop codon was most common as a stop codon. To verify the gene order of mitogenome, multiple sequence alignment of the Rana mitogenome sequences was performed. There were three kinds of gene arrangement patterns in Rana. The gene recombination mainly occurred at the control region, gene of ND5 and tRNA. The results showed that the Pattern 1 was the most common type in Rana mitogenomes, the mitochondrial gene arrangement was not related to species diversification, and several independent shifts happened in evolutionary history. Climate fluctuation and environmental change may have played an essential role in the species diversification in Rana. The molecular data obtained in this study are valuable for research on the taxonomy, population genetics, and evolution of frogs in the genera Rana. In the future, sequencing more mitogenomes from various taxonomic levels will provide useful information for better understanding the evolutionary history and will contribute to further taxonomic research within Ranidae.