#ZooKeys ZooKeys 1216: 63-82 (2024) DOI: 10.3897/zookeys.1216.131847 Research Article Insights into phylogenetic relationships and gene rearrangements: complete mitogenomes of two sympatric species in the genus Rana (Anura, Ranidae) Jingfang Li'?, Mei Xie’?, Fangpeng Zhang’, Juan Shu’, Jun Zhang’, Zhinuo Zhang?, Hongmei Xiang2®, Wansheng Jiang!2& 1 College of Biology and Environmental Sciences, Jishou University, Jishou 416000, China 2 National and Local United Engineering Laboratory of Integrative Utilization Technology of Eucommia ulmoides, Jishou University, Zhangjiajie 427000, China 3 Zhangjiajie National Forest Park, Zhangjiajie 427400, China Corresponding authors: Hongmei Xiang (hmxiang@163.com); Wansheng Jiang (jiangwschina@163.com) OPEN Qaccess Academic editor: Anthony Herrel Received: 12 July 2024 Accepted: 10 September 2024 Published: 21 October 2024 ZooBank: https://zoobank. org/09044CD6-2B7C-4AD0-8B0C- 074A798DEE93 Citation: Li J, Xie M, Zhang F, Shu J, Zhang J, Zhang Z, Xiang H, Jiang W (2024) Insights into phylogenetic relationships and gene rearrangements: complete mitogenomes of two sympatric species in the genus Rana (Anura, Ranidae). ZooKeys 1216: 63-82. https://doi.org/10.3897/ zookeys.1216.131847 Copyright: © Jingfang Li et al. This is an open access article distributed under terms of the Creative Commons Attribution License (Attribution 4.0 International - CC BY 4.0). Abstract Mitochondrial genomes (also known as mitogenomes) serve as valuable molecular markers and have found widespread applications in molecular biology and ecology. There is abundant sequence variation in vertebrate mitogenomes, and occasionally, they exhibit gene rearrangements. In this study, two Chinese endemic Rana species, Rana jiemuxiensis and Rana hanluica, were sequenced and analyzed to obtain their complete mitogenomes. The two species were sympatrically distributed in the Zhangjiajie Nation- al Forest Park, in Wulingyuan District, Zhangjiajie City, Hunan Province, China. The mitog- enome of R. jiemuxiensis was 17,506 bp, while that of R. hanluica was 17,505 bp, each comprising 13 protein-coding genes (PCGs), 22 transfer RNA genes (tRNAs), two ribo- somal RNA genes (rRNAs), and a non-coding control region (D-loop). The gene content, nucleotide composition, and evolutionary rates of each mitogenome were analyzed and compared with those of congeners. A phylogenetic analysis based on 22 mitogenomes in Rana revealed that the two sympatric species were in two different lineages, indicating that they were genetically separated to a certain extent. Three types of gene rearrange- ment patterns were identified when examining the gene orders of the 22 Rana mitoge- nomes. Most of the species shared a second and dominant gene rearrangement pattern that originated from the first ancient pattern. A “tandem duplication — multiple deletion” hypothesis was proposed to explain the evolution of these different gene rearrangement patterns. This study provided valuable data references and enhanced our understanding of the phylogenetic implications and gene rearrangements of Rana species. Key words: China, genome, mitochondrial, phylogeny, Ranine, Zhangjiajie Introduction The amphibian genus Rana Linnaeus, 1758, commonly referred to as the “wood frog” or “brown frog,” represents an early-established group typified by the Eu- ropean wood frog, Rana temporaria Linnaeus, 1758. Over the past two decades, there has been significant taxonomic restructuring within the genus Rana and its ranine counterparts. According to the Amphibian Species of the World online 63 Jingfang Li et al.: Mitogenomes of two sympatric Rana species database, the currently reclassified Rana encompasses 52 species distributed throughout temperate Eurasia into Indochina (Frost 2024). Notably, China hosts the largest number of Rana species, totaling 28, which are extensively dispersed across almost every province of the country, ranging from the north in the Hei- longjiang River basin, to the south in Hainan Island, and from the west in the Ti- betan plateau, to the east in the East China coastal plain (AmphibiaChina 2024). Numerous scientific investigations concerning Rana have focused on taxon- omy and phylogeny, as well as the discovery of new species and insights into their ecological behaviors. The first discovery of a Rana species in China took place in the Qinling Mountains during the late 19" century, officially named as Rana chensinensis David, 1875 (Boulenger 1920). However, for a long time, the Chinese wood frog R. chensinensis was often considered a synonym or sub- species of the European wood frog R. temporaria due to their morphological similarities (Liu and Hu 1961). The introduction of new phenotypic techniques and genetic tools, such as microstructure analysis, chromosome karyotyping, cytogenetics, protein analysis, isoenzymes, and DNA sequencing, has helped clarify many problematic species, identify new species, and gradually reveal phylogenetic relationships within wood frog groups (Wu 1981). For instance, Che et al. (2007) identified four well-supported clades of East Asian wood frogs through phylogenetic reconstructions using mitochondrial genes, showing sig- nificant geographic separations in the distribution patterns of some clades. Subsequently, Zhou et al. (2017) classified Rana species in East Asia into four species groups: R. japonica Boulenger, 1879, R. maoershanensis Lu, Li & Jiang, 2007, R. chensinensis and R. amurensis Boulenger, 1886. Their research em- phasized that the R. maoershanensis and R. japonica species groups formed a monophyletic clade known as the Southern Chinese wood frogs, where the R. japonica species group could be further divided into R. japonica, R. chaoc- hiaoensis Liu, 1946, and other species of the R. longicrus Stejneger, 1898 spe- cies group, each linked to distinct geographical regions. Rana hanluica Shen, Jiang & Yang, 2007 and Rana jiemuxiensis Yan, Jiang, Chen, Fang, Jin, Li, Wang, Murphy, Che & Zhang, 2011 are two species belong- ing to the R. longicrus species group (Fei et al. 2009; Wan et al. 2020), which constitutes the largest species group within the Southern Chinese Wood Frogs. Recent studies focusing on the R. Jongicrus species group have revealed the presence of several new species in China, highlighting the incomplete under- standing of this particular group of Rana (Pope and Boring 1940; Yan et al. 2011). Although initially discovered in a localized area in Hunan Province, dis- tribution records of R. jiemuxiensis and R. hanluica exhibit different patterns. Rana jiemuxiensis has a relatively restricted distribution, being found only in the surrounding areas of its type locality, the Jiemuxi National Nature Reserve in Yuanling County (Yan et al. 2011). In contrast, the distribution range of R. han- luica extends beyond its type locality in Yangmingshan, Shuangpai County (Zhu et al. 2018), encompassing neighboring provinces such as Guizhou, Guangxi, Jiangxi, and even Zhejiang Province (AmphibiaChina 2024). Notably, R. jiemux- iensis and R. hanluica coexist in the Zhangjiajie National Forest Park, represent- ing the northernmost distribution zones and likely the only coexisting area for both species to date. The reasons for the coexistence of these closely related species are multifaceted, likely stemming from their long-term evolutionary ZooKeys 1216: 63-82 (2024), DOI: 10.3897/zookeys.1216.131847 64 Jingfang Li et al.: Mitogenomes of two sympatric Rana species history. Factors contributing to their coexistence may include reproductive iso- lation, as well as potential genetic and niche differentiations. One factor fa- cilitating their coexistence could be the non-synchronous breeding seasons. Rana jiemuxiensis typically breeds from late January to mid-March, with a peak breeding period in early March (Yan et al. 2011), whereas the breeding season of R. hanluica typically occurs in October, around the “hanlu” or “Cold Dew Fes- tival” (Xia et al. 2022), one of the 24 traditional Chinese solar terms, which also inspired its species name “hanluica” (Shen et al. 2007). The genetic differentiations between the two sympatric species, R. jiemux- iensis and R. hanluica, may play a crucial role in their coexistence, yet this as- pect has received limited attention in previous studies. Furthermore, species within the Ranoidae family typically exhibit gene rearrangements in their mitog- enomes (Igawa et al. 2008), which could further contribute to genetic differen- tiation and potentially influence their coexistence. In this study, we employed next-generation sequencing technology to sequence and characterize the com- plete mitochondrial genome of R. jiemuxiensis and R. hanluica. We conduct- ed comparative analyses of these mitogenomes with those of closely related species, focusing on mitochondrial structures and features such as nucleotide composition, codon usage, and selection pressures. Additionally, we recon- structed the phylogenetic relationships using all available mitogenomes of the genus Rana obtained from NCBI and analyzed gene rearrangements within this group. The primary objective of this study is to provide novel insights into the phylogenetic implications and gene rearrangements of the two sympatric spe- cies, as well as other species within the genus Rana. Materials and methods Sample collection and sequencing Samples of R. jiemuxiensis and R. hanluica were collected from Zhangjiajie National Forest Park (29°9'39"N, 110°24'58"E), located in the Wulingyuan Dis- trict of Zhangjiajie City, Hunan Province, China. Permissions for the field survey were obtained for scientific purposes from the local administrations, and the sample collections and experimental protocols were approved by the Biomed- ical Ethics Committee of Jishou University (Approval No: JSDX-2024-0083). In accordance with the “3R principle” (Reduction, Replacement, and Refine- ment) as required by the National Ministry of Science and Technology (No. 398 [2006]), only one sample of each species was utilized. Specimens were euth- anized humanely and preserved in 85% ethanol as voucher specimens. These specimens were deposited at the Molecular Ecology Laboratory, Zhangjiajie Campus, Jishou University (R. jiemuxiensis, voucher no. JWS20211037; R. han- luica, voucher no. JWS20211131). A small volume of liver tissue was used for molecular experiments. Total DNA was extracted using the DNeasy Blood & Tissue Kit (Qiagen, Hilden, Germany). DNA library construction was performed using the VAHTS Universal DNA Library Prep Kit for Illumina V3 (Vazyme, Nan- jing, China). High-throughput sequencing was conducted on the DNBSEQ-T7 platform (Complete Genomics and MGI Tech, Shenzhen, China), generating ap- proximately 30 Gb of raw reads with a read length of 150 bp for each sample. ZooKeys 1216: 63-82 (2024), DOI: 10.3897/zookeys.1216.131847 65 Jingfang Li et al.: Mitogenomes of two sympatric Rana species Mitogenome assembly, annotation, and character analyses The complete mitogenomes of R. jiemuxiensis and R. hanluica were assem- bled using NOVOPIlasty 4.3 (Dierckxsens et al. 2017), based on the raw reads generated through next-generation sequencing (NGS). The complete mitoge- nome and the CYTB gene (1140 bp) of R. chensinensis (NCBI accession No. KF898356) were used as the “reference” and “seed” sequences for assembling the mitogenomes of both R. jiemuxiensis and R. hanluica. The positions and orientations of protein-coding genes (PCGs), ribosomal RNA genes (rRNAs), transfer RNA genes (tRNAs), and the control region (D-loop) were annotated us- ing the MITOS Web Server (http://mitos. bioinf.uni-leipzig.de/overview.py) (Ber- nt et al. 2013). The site and secondary structure of tRNAs were predicted us- ing the tRNAscan-SE 1.21 online tool (http://lowelab.ucsc.edu/tRNAscan-SE/) (Chan and Lowe 2019). Visualization and circularization of the complete mitog- enomes were performed using the CGview online server (htips://cgview.ca/) (Grant and Stothard 2008). The numbers of observed transitions (s) and transversions (v) for all the PCGs were plotted using DAMBE v7.3.11 (Xia 2018). Saturation was deter- mined based on the values of /ss (simple index of substitution saturation) and Iss.c (critical Iss value). Other analyses, such as nucleotide composition, AT and GC skewing, codon usage frequency of PCGs, and determination of se- quence genetic distances under the Kimura two-parameter (K2P) model, were conducted in MEGA 11.0. Relative synonymous codon usage (RSCU), nucle- otide diversity (Pi), and the ratio of nonsynonymous substitution rate (Ka) to synonymous substitution rate (Ks) were calculated using DnaSP 6 (Rozas et al. 2017), with all stop codons removed before analysis. Phylogenetic analysis The complete mitogenome sequences of the available species of the genus Rana were downloaded from NCBI. Twenty-two species in Rana were involved in the final dataset, including R. jiemuxiensis and R. hanluica that we sequenced. This dataset represented the most comprehensive set of mitogenome se- quences available to date. An additional species, Odorrana jingdongensis Fei, Ye & Li, 2001, was selected as the outgroup. Each of the 13 PCGs was extracted from the dataset of 23 mitogenomes and checked manually. Subsequently, all PCGs were aligned using the inbuilt MUSCLE module in MEGA and then concat- enated to create a combined PCGs dataset. The 13 PCGs concatenated data- set was used to reconstruct the phylogenetic tree using Bayesian inference (BI) (Huelsenbeck and Ronquist 2001) and maximum likelihood (ML) (H6éhler et al. 2022) methods. The optimal partition scheme and model for BI and ML were identified using PartitionFinder 2 (Lanfear et al. 2017). BI analysis was performed in MrBayes 3.2.6 with a run set to 10 million generations, sampled every 1000 generations. The initial 25% of the tree topologies were discarded as burn-in, and the consensus tree and posterior probabilities were calculated from the remaining trees. ML analysis was conducted using RAXxML 8.2.0, and the support values of the tree were assessed by conducting a bootstrap test with 1000 replicates. The phylogenetic tree was visualized, checked, and im- proved using FigTree 1.4.2 (Rambaut 2009). ZooKeys 1216: 63-82 (2024), DOI: 10.3897/zookeys.1216.131847 66 Jingfang Li et al.: Mitogenomes of two sympatric Rana species Gene rearrangement analysis Gene rearrangement in Rana was analyzed by comparing the gene orders across the entire mitogenomes of the 22 species in Rana, most of which were assembled using NGS techniques (Van Dijk et al. 2014; Ye et al. 2014). The two species we assembled in this study were compared in detail with each oth- er and within the broader context of gene rearrangement patterns in the Rana mitogenome dataset. Based on the complete mitogenomes downloaded from GenBank, we reidentified the tRNAs of each species using tRNAscan-SE (Chan and Lowe 2019) to further verify the gene orders of tRNAs. TBtools (Chen et al. 2023) was also employed to assist in the gene rearrangement analysis. Results Mitogenome assembly and annotation The complete mitogenomes of R. jiemuxiensis and R. hanluica were circular DNA molecules with lengths of 17,506 bp and 17,505 bp, respectively (Fig. 1). Both species exhibited a typical mitogenome organization, consisting of 37 genes, in- cluding 13 PCGs, 22 tRNA genes, two rRNA genes, and one control region (CR) (Table 1). All genes were encoded by the heavy strand (H-strand), except for eight tRNA genes (tRNA"®, tRNAPre, tRNAP®, tRNA‘, tRNA, tRNAS&", tRNAs, tRNA") and one PCG (ND6), which were encoded by the light strand (L-strand). The final complete mitogenomes, along with annotated information for both species, have been deposited in GenBank under accession numbers PP228843 and PP228844. In R. jiemuxiensis and R. hanluica, 13 distinct but overlapping sites were found in their mitogenomes. Three of these overlapping sites were observed between contiguous PCGs, ATP8 and ATP6, ATP6 and COX3, and NDAL and ND4, respec- tively. Ten gene intervals were also identified, with the largest one located be- tween ND5 and ND6 in both species, measuring 624 bp and 305 bp, respectively. 2 rn a jon 125 FRNA B en il ' control region _- Ae ie ' cigs ste oe ‘ ‘ ‘ a! nee _ SS = wa i Hf Wey ww hy yond Wl; May in wy Vy ty % yl cate mn ”) > hy 47; Rana hanluica 17,505bp see . Rana jiemuxiensis . 17,506bp ee ce, a e; a - 4 Phy 0 kbp 8 top 12 NW dol ie EM LH m \ hii al i PT lll sie : £1018 aie ‘ Mi trna WM rrna Meds ial d-loop | Ge Content Micc skew+ GC Skew- Zhi, ’ Wy, Teh i | Phy, iia A SY Mdiig TLL inn " Sse Figure 1. Mitodondrial gene maps of AR. jiemuxiensis and B R. hanluica. ZooKeys 1216: 63-82 (2024), DOI: 10.3897/zookeys.1216.131847 67 Jingfang Li et al.: Mitogenomes of two sympatric Rana species Table 1. Characteristics of the mitogenomes of R. jiemuxiensis (RJ) and R. hanluica (RH). Gene tRNA‘ (CUN) tRNA™ tRNA"? tRNA®re 12S rRNA tRNA“?! 16S rRNA tRNA‘ (UUR) ND1 tRNA": tRNAS&" tRNA™* ND2 tRNA"? tRNA‘? tRNA‘ NCR tRNAs tRNA™ COX1 tRNAS®" (UCN) tRNA“sP COX2 tRNA (siAsn) ATP8 ATP6 COX3 tRNACY ND3 tRNA*9 ND4L ND4 tRNAs tRNAS*' (AGY) NDS ND6 tRNAS"Y CYTB D-Loop tS 145 215 287 1217 1286 2861 2982 3895 3967 4036 4106 5139 5208 5279 5354 5377 5443 5511 7055 7128 7332 7885 7955 8110 8792 9575 9599 9984 10054 10332 11692 11760 11859 14271 14766 14838 15981 Position Length Codons i aleea RH Strand’ RJ RH lacks To Eien To RJ RH Start Stop Start Stop RJ RH codons | codons | codons | codons 74 1 74 74 74 H = = = = TAG 0 0 144 75 144 70 70 H = a = = TGT 0 0 213 145 213 69 69 L = = = 3 TGG 1 1 286 215 285 72 al H ar a = = GAA 0 0 1216 286 1215 | 930 930 H = = = = = 0 0 1285 | 1216 | 1284 69 69 H = = = = TAC 0 0 2860 | 1285 | 2859 | 1575 | 1575 H = = a = = 0 0 2934 | 2860 2933 74 74 H = = = = TAA 47 47 3935 | 2981 | 3934 | 954 | 954 H ATT AGG ATT AGG = -41 -41 3965 3894 3964 71 71 H = = = = GAT 1 0 4037 | 3965 | 4035 71 71 L = = + = TAA =2 =2 4106 4034 | 4104 71 Zul H = a = = CAT =i =i 5140 | 4104 | 5138 | 1035 | 1035 H ATG TAG ATG TAG = =2 =2 5208 | 5137 | 5206 70 70 H = = = = TCA =1 =i 5278 | 5206 | 5276 71 71 i = = > “ TGC 0 no51l 5277 | 5349 Tie 73 L = 2 = = GTT 2 2 5S28r| S352 | S374 25 26 H = = = ra = =2, =2 5442 | 5376 | 5441 66 66 L = = = = GCA 0 0 5509 | 5442 | 5508 67 67 F = = = = GTA 1 1 7064 | 5510 | 7063 | 1554 | 1554 H GTG AGG GTG AGG = coil Ad V2 8) 7126 | 7054 | 7125 72 72 L = 7 = = TGA 1 1 7196 | 7127 | 7195 69 69 H = = = = GTC 135 0 7884 | 7196 | 7883 | 553 688 H ATG T(AA) ATG T(AA) = 0 7953 | 7884 | 7952 69 69 H = = = 7 a0 a 1 1 8117 | 7954 | 8115 | 163 162 H ATG TAA ATG TAA = 8 if 8823 8109 8791 714 | 683 H ATG AGT ATG AGT ie “82 =il 9576 | 8791 | 9575 | 785 | 785 H ATG TA(A) ATG TA(A) ms =2 32 9644 | 9574 | 9643 70 70 H = = = = TCC -46 | -46 9948 | 9598 | 9947 | 350 350 H ATG TA(A) ATG TA(A) 35 35 10053 | 9983 | 10052 | 70 70 H = = = = TCG 0 0 10338 | 10053 | 10337 | 285 | 285 H ATG TAA ATG TAA = <7 <7 11703 | 10331 | 11702 | 1372 | 1372 H ATG T(AA) ATG T(AA) = =a) 12 11759 | 11691 11758 | 68 68 H 5 7 7s = GTG 0 0 11826 | 11759 | 11825 | 67 67 H = 5 = = GCT 32 31 13646 | 11857 | 13644 | 1788 | 1788 H ATG AGG ATG AGG zs 624 | 305 14765 13950 14444 | 495 | 495 L ATG AGA ATG AGA = 0 0 14834 14445 14513) 69 69 L = = = = TTC 3 3 15980 14517 15659 | 1143 | 1143 H ATG TAA ATG TAA = 0 0 17505 | 15660 | 17505 | 1525 | 1846 H az 5 = = = 0 0 * H and L indicate genes transcribed on the heavy and light strand, respectively. # Positive numbers correspond to the nucleotides separating adjacent genes; negative numbers indicate overlapping nucleotides. The shortest interval identified was only 1 bp, found in multiple locations within both species. The lengths of the 13 PCGs varied considerably. The longest gene was ND5, which was identical in length in both species at 1788 bp, while the shortest gene was ATP8, with lengths of 162 bp and 163 bp for R. jiemuxiensis and R. hanluica, respectively. Both species had ATG as the start codon for most PCGs, except for ND1 (ATT) and COX1 (GTG). Five typical stop codons, includ- ing TAG, AGG, AGA, AGT, and TAA, as well as two kinds of incomplete terminal codons (TA-, T-), were found in the PCGs within their mitogenomes. ZooKeys 1216: 63-82 (2024), DOI: 10.3897/zookeys.1216.131847 68 Jingfang Li et al.: Mitogenomes of two sympatric Rana species Nucleotide composition and diversity The overall base composition of R. jiemuxiensis was as follows: A (24.33%), T (29.11%), G (15.65%), C (30.49%), while that of R. hanluica was A (24.47%), T (29.37%), G (15.65%), C (30.49%). Both species showed an A+T bias with great- er A+T than G+C content. Additionally, both species exhibited negative AT skew and GC skew, indicating a predominant bias towards T and C base pairs. The mitogenome sequences of the 22 Rana species compiled in this study ranged from approximately 16,000 bp to 22,000 bp in length, indicating a complex mi- togenome evolution among Rana species. However, all species showed a simi- lar A+T content bias and T and C base pair biases that resembled R. jiemuxien- sis and R. hanluica (Table 2). When examining the average length and nucleotide composition of each PCG (Table 3), there were generally similarities, but differences remained. The shortest PCG was ATP8, and the longest one was ND5. The ND6 gene exhib- ited distinct differences in AT skew and GC skew compared to other PCGs. After removing the stop codon from each PCG, the total aligned length of the final 13 PCGs dataset was 11,244 bp. Nucleotide diversity analysis of each PCG showed values ranging from 0.18 (COX3) to 0.27 (ND5). In addition to ND5 (0.27), four other PCGs exhibited relatively high nucleotide diversity, name- ly ND3 (0.24), ATP6 (0.24), ND2 (0.25), and ATP8 (0.26), while the remaining PCGs showed relatively low nucleotide diversity, all less than 0.2. Table 2. Basal composition (percentage) of the mitogenomes of R. jiemuxiensis and R. hanluica and 20 other Rana species. R. R. R. R. R. R. R R. R. R. DD VY DVD VDD VD D VD Name hanluica jiemuxiensis dybowskii chensinensis draytonii huanrensis . amurensis chaochiaoensis kukunoris omeimontis temporaria uenoi johnsi dabieshanensis hanluica zhenhaiensis wuyiensis catesbeiana arvalis coreana longicrus kunyuensis Avg. *, the sequence obtained in this study. ZooKeys 1216: 63-82 (2024), DOI: 10.3897/zookeys.1216.131847 T% 29.37567 29:19775 30.53428 30.69457 29.67805 30.76512 30.9651 30.19221 30.70901 29.75714 30.66239 30.60552 29.40915 29.61744 29.37461 29.16111 29.43584 32.8264 30.69122 29.22448 31.33066 31.27726 30.24552 C% A% 6% Total length — (A+T)% GC skew | ATskew s edsca 30.49626 | 24.47528 | 15.65279 17505 53.85094 | -0.32164 | -0.091 | PP228844% 30.73639 | 24.33298 | 15.81288 17506 53.45073 | -0.3206 | -0.08952 | PP228843* 29.31434 | 2457703 | 15.57435 18864 55.11131 | -0.30609 | -0.1081 | KF898355 29.10062 | 24.94212 15.26269 18808 55.63669 | -0.31192 | -0.10339 | KF898356 30.06937 | 25.37353 | 14.87905 17805 55.05158 | -0.33795 | -0.07819 | KP013110 29.04605 | 25.06458 | 15.12425 19253 55.8297 | -0.31518 | -0.10211 | KT588071 29.11325 | 25.5609 | 1436075 18470 56.526 | -0.33934 | -0.09561 | KU343216 29.76508 | 2468411 | 15.3586 18591 54.87631 | -0.31927 | -0.10037 | KU246048 29.11663 | 25.06005 | 15.11431 18863 55.76906 | -0.31657 | -0.10129 | KU246049 30.03292 | 2416155 | 16,04839 19934 53.91869 | -0.30347 | -0.10378 | KU246050 29.20228 | 24.90207 15.23326 16061 55.56446 | -0.31437 | -0.10367 | MH536744 29.439 24.87979 | 15.07569 17370 55.48531 | -0.32266 | -0.10319 | Mwo09067 30.72611 | 25.2981 | 14.56665 17837 54.70724 | -0.35678 | -0.07515 | MZ571365 30.22242 24.1637 ~—-15.99644 18291 53.78114 | -0.3078 | -0.10141 | Mw526989 30.5133 24.4907 | 15.62139 19395 53.86531 | -0.32279 | -0.09067 | MZ680529 30.59336 | 24.66862 | 15.57691 18806 53.82973 | -0.32524 | -0.08346 | OL681880 30.66382 | 25.44937 | 14.45097 17779 54.88521 | -0.35937 | -0.07263 | 0L467321 26.68979 | 25.99609 | 14.48773 17212 58.82248 | -0.29633 | -0.11612 | ON746668 29.13442 | 24.81986 15.35451 16143 55.51108 | -0.30974 | -0.10577 | MT872666 30.46958 | 24.7243 | 15.58164 22262 53.94877 | -0.32329 | -0.08342 | O0N920705 28.63373 | 25.73209 14.30352 17833 57.06275 | -0.33375 | -0.09811 | MZ680528 28.63373 | 25.8834 | 14.20561 22255 57.16066 | -0.3368 | -0.09436 | KF840516 29.62343 | 24.96542 | 15.16564 18493 55.21094 | -0.3228 | -0.09564 2 69 Jingfang Li et al.: Mitogenomes of two sympatric Rana species Table 3. The average proportion of 13 PCGs in 22 Rana species. Gene ND1 ND2 COX1 COX2 ATP8 ATP6 COX3 ND3 ND4L ND4 ND5 ND6 CYTB T% 31.8 29 29.7 26.3 28.9 31.5 30.1 33:2 30.1 30:3 30.2 34.8 29.2 C% 31.5 32 28 27.5 28.8 SA. 30.4 31.3 31.8 31 29.8 10.7 32.8 A% G% (A+T)% Total length GC skew AT skew 23.7 13 55.5 958 -0.41964 -0.14595 27.7 11.3 56.7 1029 -0.43921 -0.02293 24.8 17.4 54.5 1684 -0.26115 -0.08991 29.7 16.5 56 548 -0.22897 0.060714 32.1 10.1 61 159 -0.48205 0.052459 25.4 11.2 56.9 682 -0.47541 -0.10721 22.7 16.8 52.8 783 -0.28358 -0.14015 21 14.6 54.2 337 -0.38912 -0.22509 24.6 13.5 54.7 282 -0.38073 -0.10055 25.9 12.9 56.2 1363 -0.40278 -0.07829 25.7 14.3 55.9 1779 -0.3573 -0.0805 17.8 36.6 52.6 491 0.02521 -0.32319 24.1 13.9 53.3 1140 -0.35499 -0.09568 Analysis of codon usage and genetic distance There are a total of 20 amino acids encoded by the PCGs of R. jiemuxiensis and R. hanluica. Among these amino acids, Leu, Ser, and Arg had the highest frequency, while Trp and Met had the lowest. According to the RSCU analysis, Leu, Ser, and Arg were encoded by six codons each; Pro, Thr, Val, Ala, and Gly were encoded by four codons each; Phe, Tyr, Cys, His, Gln, Asn, Lys, Asp, and Glu were encoded by two codons each; Trp and Met were encoded by only one codon each. This reflects a significant bias in codon usage in their mitoge- nomes (Fig. 2). Genetic distances were analyzed among 22 Rana species (Fig. 3). Rana jiemuxiensis was found to be closely related to R. zhenhaiensis Ye, Fei & Mat- sui, 1995, R. longicrus, R. hanluica, R. dabieshanensis Wang, Qian, Zhang, Guo, Pan, Wu, Wang & Zhang, 2017, and R. omeimontis Ye & Fei, 1993, with genetic distances of 0.08327, 0.08529, 0.08547, 0.08932, and 0.09243, respectively. Conversely, R. hanluica was closely related to R. dabieshanensis, R. omeimon- tis, R. jiemuxiensis, R. zhenhaiensis, and R. longicrus, with genetic distances of 0.07802, 0.07793, 0.08357, 0.08398, and 0.08609, respectively. Rana catesbe- iana Dubois, 1992 appeared to be the most genetically distant from all other Rana species, suggesting it may be an ancient ancestor species. Sliding win- dow analysis along the concatenated PCGs dataset showed significant varia- tion in nucleotide diversity (Pi) among different genes (Fig. 4A). Substitution saturation test indicated that /ss < /ss.c in general, and scatter diagrams also suggested that the PCGs substitution was not saturated, making them suitable for phylogenetic analysis (Fig. 4B). Standard deviations of Ka and Ks for the 13 PCGs across 22 species showed that the data were generally concentrated, with variances of Ks generally great- er than those of Ka (Fig. 5A). ATP8 and ND6 exhibited the highest and low- est standard deviations of Ks, respectively, while ND5 and COX1 showed the highest and lowest standard deviations of Ka, respectively. The average Ka/Ks values from pairs of species among the 22 Rana species fell into the range of 0 to 1. Among the 13 PCGs analyzed, ND6 and ATP8 genes evolved relatively quickly, exhibiting the highest Ka/Ks values. Conversely, COX1 and CYTB had the lowest Ka/Ks ratios (Fig. 5B). However, all 13 PCGs had Ka/Ks values below 0.45, with no signals of positive selection detected. ZooKeys 1216: 63-82 (2024), DOI: 10.3897/zookeys.1216.131847 70 Jingfang Li et al.: Mitogenomes of two sympatric Rana species o> [o> oa o & & w ioe) ho pe) = Relitive Synonymous Codon Usage (RSCU) Peace: Synonymous Codon Usage (RSCU) Oo Phe Leu Ser Tyr Cys Trp Pro His Gln Arg Ile Met Thr Asn Lys Val Ala Asp Glu Gly Phe Leu Ser Tyr Cys Trp Pro His Gln Arg Ile Met Thr Asn Lys Val Ala Asp Glu Gly Figure 2. Relative synonymous codon usage in the protein coding genes of A R. jiemuxiensis and B R. hanluica. R. hanluica e R. jiemuxiensis R. dybowskii R. chensinensis R. draytonii of 0.15- R. huanrensis 0.1- R. amurensis R. chaochiaoensis alia 0.05 R. kukunoris a R. omeimontis R. temporaria R. uenoi R. johnsi R. hanluica R. zhenhaiensis R. wuyiensis R. catesbeiana R. arvalis R. longicrus R. kunyuensis R. coreana 123456789 10 11 12 13 14 15 16 17 18 19 20 21 Figure 3. Genetic distance heatmap plot within 22 Rana species. Phylogenetic analysis The tree topologies resulting from BI and ML analyses were identical, with only slight differences in the support values of some nodes (Fig. 6). Generally, the support values were high in most branches. Despite their sympatric coexistence, R. jiemuxiensis and R. hanluica did not cluster together in the same subclade. Rana jiemuxiensis was grouped with R. longicrus and R. zhenhaiensis to form one subclade, while R. hanluica was grouped with R. omeimontis and R. dabieshanen- sis in another subclade. However, the subclades containing R. jiemuxiensis and ZooKeys 1216: 63-82 (2024), DOI: 10.3897/zookeys.1216.131847 71 Jingfang Li et al.: Mitogenomes of two sympatric Rana species 4 tw oN ’ Ay : he 0. 03 s and v Av 0. 00 0.0000 0.0461 0.0921 0.1382 0.1842 0.2303 0.2764 Figure 4. A Nucleotide diversities and B substitution saturation plot of the mitochondrial protein coding genes. Ka/Ks Figure 5. Standard deviation of A Ks and Ka and B the Ka/Ks ratio of 13 protein coding genes. R. hanluica were then grouped together to form a major clade. The phylogenetic tree also revealed five other major clades, with R. catesbeiana and R. draytonii Baird & Girard, 1852 representing the two most ancient clades. Evolutionary branch lengths reflected the evolutionary history of each branch, indicating that earlier clades had longer branch lengths. Sequence lengths, distribution of alti- tude, and degree of threatened levels, however, did not show strong links to the phylogenetic relationships, indicating a complicated and specific evolutionary history for each species. Gene rearrangement analysis As expected, the mitogenomes of Rana species exhibited substantial gene rear- rangements and were categorized into three distinct patterns (Fig. 7A). Only two species, R. draytonii and R. zhenhaiensis, retained the first pattern, which is as- sumed to be ancient based on its presence in most amphibian species, including the relatively close outgroup species Odorrana jingdongensis used in this study, as well as species from other studies such as Leptobrachium and Boulenophrys species in Anura (Zhou et al. 2023; Xiang et al. 2024) and Tylototriton species ZooKeys 1216: 63-82 (2024), DOI: 10.3897/zookeys.1216.131847 72 Jingfang Li et al.: Mitogenomes of two sympatric Rana species Degree of threats(DT) @ ww O LC ® DD =) NE Distribution of altitude(DA) | Over 2000meters | 1000-2000meters ia 500-1000meters [] Under 500meters Tree scale: 0.1. —— bootstrap B vs) - B < oO 4 oO > ws-+-+--5- -Rana omeimontis (KU246050) ore 19934 O HA si ~-----~=~=-Rana dabieshanensis (MW526989) yous aids O TM eg wn ee ee Rana hanluica (MZ680529) Ponor Zone O a 9 0lCl (tttt*«*«i‘“ ‘é*‘CS Rana hantutca (PP228844) A 17505 @® H crecece Rana longicrus ¢ MZ680528 ) 0.03 raid @ an OD ee ‘Rana thenhaiensis (OL681880) 0.039 10S O L] 1000 Lng 2 = = = 2 Rana jiemuxiensis (PP228843) men 17505 'S) E] ~ = === ===» ~~ -Rana chaochiaoensis (KU246048) i 0.188 a @ [O ——— 0.357 a @ a] aes 18808 O | i 0.065 x ) an _— 18864 @ oO a 0.213 on @ | a 17833 O Be ae 22255 oO El a, , 18470 @ an ~ anes esse sss sarass ‘Rana temporaria (MH536744 ) aoe ioe O B wart tt SoS eR SSS ee Rane ama ery pcooe i 0210 els O [] is - 7s uae Rana johnsi (MZ571365 ) =a e e > - 22 Seen -Rana wuyiensis (OL467321) — @ a lubiatticoe 8 Rana draytonii_ (KPO013110) =a e ee — O (ix) e —iiniiiitinir ini ein = 5 5 = -Odorrana jingdongensis (MT757792) sc 0.812 Figure 6. Phylogenetic tree of Rana species based on BI and ML analyses. Note: samples sequenced in this study are high- lighted in red. BRL represents the length of the evolutionary branch and NL represents the length of the nucleotide sequence. in Caudata (Wang et al. 2022). Interestingly, 18 out of the 22 examined Rana species shared the second and also the most dominant pattern, accounting for 82% of the species analyzed. Both R. jiemuxiensis and R. hanluica, sequenced in this study, belonged to this dominant pattern. Two other species, R. amuren- sis Boulenger, 1886 and R. coreana Okada, 1928, shared the third, rare pattern, which is similar to the second pattern except that the ND5 gene is transposed to a position behind the control region (D-loop). The transition from the first pattern to the second and third patterns primarily resulted from changes in the positions of tRNA genes. The gene orders of rRNAs remained unchanged, except for the rearrangement of ND5 from the second to the third pattern. We proposed a plausible scenario to explain the mitochondrial gene rear- rangements within Rana, considering that duplications and losses are more like- ly to occur among tRNAs than rRNAs and PCGs. According to the principle of parsimony, a “tandem duplication - multiple deletion” event in a sequence region spanning from ND4 to the D-loop likely triggered the transition from the first an- cient pattern to the second pattern. Subsequently, a simple transposition of ND5 would have driven the second pattern to evolve into the third pattern (Fig. 7B). In the evolutionary history of Rana, for reasons yet unclear, the second pattern became the dominant pattern within this group. This dominance pattern in mito- chondrial gene orders in Rana is distinct from those in other amphibian species. Species verification based on individual genes Species verification is crucial for publishing the complete mitogenome of a spe- cies. To verify the molecular identification of the species involved in this study, we selected the individual 16S rRNA, ND2, and CYTB genes as target genes because ZooKeys 1216: 63-82 (2024), DOI: 10.3897/zookeys.1216.131847 73 Jingfang Li et al.: Mitogenomes of two sympatric Rana species R. omeimontis R. dabieshanensis R. hantuica ait R. hantuica R. longicrus Rzhenhaiensis — GOBOCBOCDOOOIOOOOOOHOOSOSELIGODOLDD9OCSBOBOO™ (1) R. jiemuxiensis R. chaochiaoensis R. kukunoris R. chinesensis (tl) R. huanrensis R. dybowshii R. uenoi R. coreana » (all) R, kunyuensis a R. amurensis SEO OCOROBOENSUTSSOOCOMOCMOMEBCHOEROCEOS + (IID) R. temporaria arias an R. johnsi R. wiyiensis R. draytonii BOBODQINAIOOGOOOOOOODODSOSGSOSODOISOSCBBOBOO + (1 R. catesbeiana OOOMEOSBOBOOODQOOOOCSSOOS2SDEOSCEBDCSBOB— 1) O. jingdongensis — TBOCBODOOCHOOOQOOCHOOBOSHSODOGBOSCSBOSO® + @ D-loop @ 12S@ 16s @ ND! O nv2 B cox: Bcoxn Wares Carrs B coxm G no3s Gnosa_ G nos W os @ ndoB cyts B | tandem duplication BoccEmosoraa Gooceseosoos x xx ¥ x XXX VV xX Figure 7. A Phylogenetic gene orders within 22 Rana species and B three patterns of mitochondrial gene rearrangement. Note: The icons represent tRNAs as: (1) L1: tRNA‘* (CUN), (2) T: tRNA™’, (3) P: tRNAP®, (4) F: tRNAPhe, (5) V: tRNAY2!, (6) L2: tRNA‘® (UUR), (7) I: tRNA’, (8) Q: tRNAS", (9) M: tRNA™*, (10) W: tRNA™®, (11) A: tRNA‘®, (12) N: tRNA*S", (13) O: NCR, (14) C: tRNA“, (15) Y: tRNA™, (16) S1: tRNA‘ (UCN), (17) D: tRNA**?, (18) K: tRNA“ 45»), (19) G: tRNA®", (20) R: tRNA“, (21) H: tRNAs, (22) S2: tRNAS* (AGY), (23) E: tRNAS". they have abundant resources in NCBI based on previous studies (Djebbi et al. 2018). These genes are also frequently used in DNA barcoding studies (Formen- ti et al. 2021; Ahmed et al. 2022). Through BLAST homology searching opera- tions, the CYTB, ND2, and 16S rRNA fragments from our samples showed high similarities, 100%, 99.64%, and 99.79%, respectively, with the species named R. jiemuxiensis in NCBI, which was collected from its type locality in Yuanling Coun- ty, Hunan Province. Similarly, the gene fragments of R. hanluica samples showed nearly 100% similarity with these genes of R. hanluica collected from Lishui City, Zhejiang Province. It is generally accepted that gene similarities between individ- uals of the same species are usually greater than 98%. Therefore, we conclude that the two species used in this study were correctly identified. Discussion Characteristics of the mitogenomes Mitochondrial DNAs, or mitogenomes, often serve as valuable molecular mark- ers and have been widely applied in molecular biology and ecological studies. ZooKeys 1216: 63-82 (2024), DOI: 10.3897/zookeys.1216.131847 74 Jingfang Li et al.: Mitogenomes of two sympatric Rana species Typically, animal mitogenomes contain 2 rRNAs (12S and 16S rRNA), 22 tRNAs, 13 PCGs, and a control region (also known as the D-loop), with a sequence length usually ranging from 16 to 17 kb (Huang et al. 2016; Tang et al. 2021). Vertebrate mitogenomes are particularly conservative with several unique characteristics, including maternal inheritance, a rapid evolutionary rate, and low levels of recombination. These characteristics make mitochondrial DNA valuable in reconstructing phylogenetic relationships, testing selective pres- sures, and identifying species using mitochondrial barcoding genes, etc. (Jiang et al. 2023; Lan et al. 2023; Zhang et al. 2023; Xiang et al. 2024). The utilization of mitochondrial DNA in molecular identification provides a valuable tool in taxonomic studies compared to traditional morphological approaches (Li et al. 2021). In recent years, an increasing number of species identifications have relied on the combination of morphology and molecular evidence. In some cases, such as cryptic species identifications, molecular ev- idence carries more weight than morphological evidence. Molecular data are commonly used in species identification and phylogenetic studies (Miya and Nishida 2015; Wang et al. 2016; Tan et al. 2018), especially with the rapid ad- vancements in NGS technology. The two sympatric Rana species involved in this study, R. jiemuxiensis and R. hanluica, are morphologically similar, with a distinguishing characteristic primarily based on the number of bands on the thighs and tibia (5 or 6 in R. jiemuxiensis vs 8 or 9 in R. hanluica). In the field, these two species often seem to completely overlap, as they are sometimes found coexisting in very small areas. However, their genetic distance and phy- logenetic position are distinct (Fig. 3 and Fig. 6, respectively), which is also evident in the species verification based on BLAST searches of our sequences against the NCBI database. The nucleotide composition of both R. jiemuxiensis and R. hanluica exhibited a distinct A+T rich pattern (Table 3), which is commonly observed in the mitoge- nomes of other species (Rayko 1997; Francoso et al. 2023). This typical nucleo- tide composition is highly conserved among vertebrates (Zhou et al. 2006; Hao et al. 2016; Vinogradov and Anatskaya 2017). In Rana species, the A+T content gen- erally exceeds 50% but is less than 60%. The gene structure of both sequenced species was identical and similar to other animal species, with the majority of genes, especially the PCGs, located on the H strand, while only a few are on the L strand (Boore 1999). There were 13 gene overlaps in the mitogenomes of R. jiemuxiensis and R. hanluica, which is also common in other Rana species (Zyla et al. 2019), indicating the compact nature of mitogenomes in regulating gene expression. However, gene intervals were also commonly present in frog mitog- enomes, with the locations varying more than overlapping regions and showing species-specific patterns. For example, for the two Rana species here, the gene interval between tRNA" and ND1 was 41 bp, whereas there was no interval be- tween 12S rRNA and tRNA‘? In contrast, in the mitogenomes of two Leptobrachi- um species (Zhou et al. 2023), there was no interval between tRNA‘ and ND1 genes, and the interval between 12S rRNA and tRNA“ was 4 bp. The usage of start and stop codons in different mitogenomes also showed species-specific patterns, which have received limited attention (Wang et al. 2022). In the two Rana mitogenomes sequenced here, all genes shared the start codon ATG, ex- cept for COX1 (GTG) and ND1 (ATT). In comparison, the start codons of COX1, ATP6, and ATP8 were GTG in the mitogenomes of two Tylototriton species. Three ZooKeys 1216: 63-82 (2024), DOI: 10.3897/zookeys.1216.131847 75 Jingfang Li et al.: Mitogenomes of two sympatric Rana species rare types of stop codons (AGT, AGA, and AGG) were also identified in this study, which are less common but have been observed in other Rana species. The ratio of Ka and Ks is a popular proxy for detecting adaptive evolution, with Ka/Ks > 1 reported in the mitochondrial PCGs of some species (Ye et al. 2017; Bi et al. 2020). However, none of the Ka/Ks values from the PCGs within Rana species exceeded 1 (Fig. 5A), indicating that the overall evolution pattern of Rana mitogenomes tends to be conservative in maintaining the functions of regularly generated proteins. The relatively high Ka/Ks ratios observed in some PCGs may represent higher evolutionary rates, such as those of the ATP6, ND3, ND5, and ND6 genes (Fig. 5B). These mitochondrial PCGs that evolve more rapidly may accumulate advantageous mutations, potentially enhancing the fitness of the species in adapting to environmental changes (Yang et al. 2018). Phylogenetic implications and gene rearrangement evolution Mitochondrial DNA sequences, particularly mitogenomes, are increasingly utilized in phylogenetic studies (Dhorne-Pollet et al. 2020). In this study, a phylogenetic tree of the genus Rana was constructed based on 13 mitochondrial PCGs from 22 species, representing the most comprehensive and detailed mitogenomic phy- logenetic tree of Rana to date. Although coexisting in Zhangjiajie National Forest Park, R. jiemuxiensis and R. hanluica were found in distinct phylogenetic lineages (Fig. 6). Rana jiemuxiensis clustered with R. longicrus and R. zhenhaiensis to form a subclade, while R. hanluica grouped with R. omeimontis and R. dabieshanensis in another subclade. This supports our speculation that the two sympatric species are genetically separated to a certain extent. Despite including more species, our phylogenetic tree was largely consistent with previous studies that also utilized mitogenomes in Rana (Chen et al. 2018; Wang et al. 2020; Zhang et al. 2020). The differentiation of species in Rana is possibly associated with geographic isolation and habitat selection. The known distribution altitudes of Rana spe- cies were mapped onto the phylogenetic tree (Fig. 6). It revealed that Rana spe- cies have a wide altitudinal adaptation, ranging from over 100 meters to more than 2000 meters above sea level. However, no clear correlations were evident between phylogenetic position and distribution altitude ranges. Previous stud- ies have shown that the skin morphological properties of Rana species at high altitudes differ from those at low altitudes (Zhi and Li 2016). It was further inferred that amphibians living at high altitudes, thriving in low temperatures, may have slower evaporative rates (Zyla et al. 2019). Additionally, the degree of threat for each Rana species, acquired from the IUCN Red List of Threatened Species of Amphibians, was also mapped onto the phylogenetic tree. Among the 22 Rana species in this study, five were classified as vulnerable (VN), 15 as least concern (LC), three as data deficient (DD), and four as not evaluated (NE). Generally, the species in Rana are experiencing a medium level of threats. Previous studies have revealed that interspecies variations in mitochondrial gene orders are prone to occur in certain groups (Shao and Barker 2003; Liu et al. 2017), including species within Rana (Zhu et al. 2018). However, few studies have analyzed mitogenomic gene rearrangements in Rana from a phylogenetic perspective. By examining the gene orders of 22 Rana species in this study, three distinct patterns can be discerned (Fig. 7A). The first pattern represents a very ancient type that is shared with most amphibian species (Wang et al. 2022; Zhou ZooKeys 1216: 63-82 (2024), DOI: 10.3897/zookeys.1216.131847 76 Jingfang Li et al.: Mitogenomes of two sympatric Rana species et al. 2023). The second pattern, however, is dominant in Rana species and mainly results from position changes of some tRNAs in their mitogenomes. The third pat- tern appears to be derived from the second one, with only one PCG (ND5) seem- ingly transposed. Interestingly, the three patterns of gene arrangements show a trend of parallel evolution, with none of them being confined to a single subclade. A simple parsimony scenario is proposed to explain the evolutionary process of the three patterns observed in this study (Fig. 7B), which we refer to as the “tandem duplication-multiple deletion” hypothesis. This phenomenon is com- monly observed in mitogenomes, where tandem duplications can occur due to strand slippage during replication, followed by gene deletions under selective pressures (Moritz and Brown 1987; Levinson and Gutman 1987). In the case of Rana mitochondrial genes, this duplication and loss event can drive the evolu- tion from the first ancient pattern to the second dominant pattern. The genes involved in this process include 6 tRNAs, 4 PCGs, and a control region, with changes mainly resulting from the tRNAs. Previous studies have indicated that tRNAs in mitogenomes are relatively neutral (Dowton et al. 2009; Lavrov and Pett 2016), suggesting that gene arrangements resulting from tRNA position changes may not significantly alter the conservative functions of mitochondria. Additionally, tandemly duplicated regions often favor genes that utilize stem- loop structures during replication (Stanton et al. 1994), a characteristic typical of vertebrate tRNAs (Macey et al. 1997). The involvement of one PCG (ND5) from the second to the third pattern echoes that changes in tRNAs and PCG po- sitions are frequently observed in mitochondrial gene rearrangements in both invertebrates and vertebrates (Cantatore et al. 1989; Miya and Nishida 1999). The study of gene rearrangement patterns can provide insights into taxon- omy and elucidate the complex evolutionary history among species (Ye et al. 2021; Feng et al. 2023). However, a comprehensive phylogenetic tree compris- ing most species is necessary to fully understand the evolutionary background that underlies the gene rearrangements in Rana. As more mitogenomes, such as those of the two species in this study, become available in the future, we antic- ipate that the complete evolutionary history of Rana will gradually be unveiled. Additional information Conflict of interest The authors have declared that no competing interests exist. Ethical statement The collection and handling of Rana species in this study were approved by the Biomed- ical Ethics Committee of Jishou University (Approval No: JSDX-2024-0083). Funding This work was supported by the Graduate Scientific Research Innovation Project of Hu- nan Province (CX20231083), National Natural Science Foundation of China (32060128), and Zhilan Foundation (2020040371B/20220100118B). Author contributions Conceptualization, H.X. and W.J.; methodology, J.L. and W.J.; software, M.X. and F.Z.; validation, W.J. and J.L.; formal analysis, J.L. and W.J.; investigation, M.X., F.Z. J.S., J.Z. ZooKeys 1216: 63-82 (2024), DOI: 10.3897/zookeys.1216.131847 77 Jingfang Li et al.: Mitogenomes of two sympatric Rana species and Z.Z.; resources, J.S., J.Z. and Z.Z.; data curation, J.L.; writing original draft prepara- tion, J.L.; writing review and editing, H.X. and W.J.; visualization, J.L.; supervision, W.J.; project administration, W.J.; funding acquisition, W.J. All authors have read and agreed to the published version of the manuscript. Author ORCIDs Hongmei Xiang ® https://orcid.org/0009-0006-3251-3726 Wansheng Jiang © https://orcid.org/0000-0002-6498-944X Data availability The final complete mitogenomes, along with annotated information for both species, have been deposited in GenBank under accession numbers PP228843 and PP228844. All the analyses and findings of this study were based on these two sequences and oth- er sequences that were available in GenBank. References Ahmed §S, Ibrahim M, Nantasenamat C, Nisar MF, Malik AA, Waheed R, Anmed MZ, Ojha SC, Alam MK (2022) Pragmatic applications and universality of DNA barcoding for substantial organisms at species level: a review to explore a way forward. BioMed Research International 2022(1): 1846485. https://doi.org/10.1155/2022/1846485 AmphibiaChina (2024) The database of Chinese amphibians. Kunming Institute of Zool- ogy (CAS), Kunming, Yunnan, China. https://www.amphibiachina.org/ [accessed on 6 Jun 2024] Bernt M, Donath A, Juhling F, Externbrink F, Florentz C, Fritzsch G, Pitz J, Middendorf M, Stadler PF (2013) MITOS: improved de novo metazoan mitogenome annotation. Molecular Phy- logenetics and Evolution 69(2): 313-319. https://doi.org/10.1016/j.ympev.2012.08.023 BiC, LUN, Xu Y, He C, Lu Z (2020) Characterization and analysis of the mitochondrial genome of common bean (Phaseolus vulgaris) by comparative genomic approaches. Internation- al Journal of Molecular Sciences 21(11): 3778. https://doi.org/10.3390/ijms21113778 Boore JL (1999) Animal mitochondrial genomes. Nucleic Acids Research 27(8): 1767- 1780. https://doi.org/10.1093/nar/27.8.1767 Boulenger GA (1920) A monograph of the South Asian, Papuan, Melanesian and Austra- lian frogs of the genus Rana. Records of the Zoological Survey of India 20(1): 1-223. https://doi.org/10.2651 5/rzsi/v20/i1/1920/163533 Cantatore P, Roberti M, Rainaldi G, Gadaleta MN, Saccone C (1989) The complete nucle- otide sequence, gene organization, and genetic code of the mitogenome of Paracen- trotus lividus. The Journal of Biological Chemistry 264(19): 10965-10975. https:// doi.org/10.1016/S0021-9258(18)60413-2 Chan PP, Lowe TM (2019) tRNAscan-SE: searching for tRNA genes in genomic se- quences. Methods in Molecular Biology (Clifton, N.J.) 1962: 1-14. https://doi. org/10.1007/978-1-4939-91 73-0_1 Che J, Pang J, Zhao EM, Matsui M, Zhang YP (2007) Phylogenetic relationships of the Chi- nese brown frogs (genus Rana) inferred from partial mitochondrial 12S and 16S rRNA gene sequences. Zoological Science 24(1): 71-80. https://doi.org/10.2108/zsj.24.71 Chen J, Xing Y, Yao W, Zhang C, Zhang Z, Jiang G, Ding Z (2018) Characterization of four new mitogenomes from Ocypodoidea & Grapsoidea, and phylomitogenomic in- sights into thoracotreme evolution. Gene 675: 27-35, 47. https://doi.org/10.1016/j. gene.2018.06.088 ZooKeys 1216: 63-82 (2024), DOI: 10.3897/zookeys.1216.131847 78 Jingfang Li et al.: Mitogenomes of two sympatric Rana species Chen C, Wu Y, Li J, Wang X, Zeng Z Xu J, Liu Y, Feng J, Chen H, He Y, Xia R (2023) TBtools- Il: a “one for all, all for one” bioinformatics platform for biological big-data mining. Molecular Plant 16(11): 1733-1742. https://doi.org/10.1016/j.molp.2023.09.010 Dhorne-Pollet S, Barrey E, Pollet N (2020) A new method for long-read sequencing of an- imal mitochondrial genomes: application to the identification of equine mitochondrial DNA variants. BMC Genomics 21(1): 785. https://doi.org/10.1186/s12864-020-07183-9 Dierckxsens N, Mardulyn P, Smits G (2017) NOVOPlasty: de novo assembly of organelle genomes from whole genome data. Nucleic Acids Research 45(4): e18. https://doi. org/10.1093/nar/gkw955 Djebbi S, Mezghani-Khemakhem M, Bouktila D, Makni H, Makni M (2018) Assessment of pea weevil Bruchus pisorum (Coleoptera, Bruchidae) genetic diversity based on mitochondrial COI gene sequences. African Entomology 26: 95-103. https://doi. org/10.4001/003.026.0095 Dowton M, Cameron SL, Dowavic JI, Austin AD, Whiting MF (2009) Characterization of 67 mitochondrial tRNA gene rearrangements in the Hymenoptera suggests that mito- chondrial tRNA gene position is selectively neutral. Molecular Biology and Evolution 26(7): 1607-1617. https://doi.org/10.1093/molbev/msp072 Fei L, Hu S, Ye C, Huang Y (2009) Anura Ranidae. Fauna Sinica; Science Press: Beijing, Chinas= 1422. Feng H, Lv S, Li R, Shi J, Wang J (2023) Mitochondrial genome comparison reveals the evolution of cnidarians. Ecology and Evolution 13(6): e10157. https://doi. org/10.22541/au.166687211.16942420/v1 Formenti G, Rhie A, Balacco J, Haase B, Mountcastle J et al. (2021) Complete vertebrate mitogenomes reveal widespread repeats and gene duplications. Genome Biology 22(1): 120. https://doi.org/10.1186/s13059-021-02336-9 Francoso E, Zuntini AR, Ricardo PC, Santos PKF, Souza Araujo N, Silva JPN, Gongalves LT, Brito R, Gloag R, Taylor BA, Harpur BA, Oldroyd B P Brown MJF, Arias MC (2023) Rapid evolution, rearrangements and whole mitogenome duplication in the Australian stingless bees Tetragonula (Hymenoptera, Apidae): a steppingstone towards under- standing mitochondrial function and evolution. International Journal of Biological Macromolecules 242(1): 124568. https://doi.org/10.1016/j.ijbiomac.2023.124568 Frost DR (2024) Amphibian Species of the World, an Online Reference. https://amphibi- ansoftheworld.amnh.org/ [accessed on 6 Jun 2024] Grant JR, Stothard P (2008) The CGView Server: a comparative genomics tool for cir- cular genomes. Nucleic Acids Research 36: W181-W184. https://doi.org/10.1093/ nar/gkn179 Hao S, Ping J, Zhang Y (2016) Complete mitochondrial genome of Gekko chinensis (Squamata, Gekkonidae). Mitochondrial DNA Part A 27(6): 4226-4227. https://doi. org/10.3109/19401 736.2015.1022751 Hohler D, Pfeiffer W, loannidis V, Stockinger H, Stamatakis A (2022) RAxML Grove: an empirical phylogenetic tree database. Bioinformatics 38(6): 1741-1742. https://doi. org/10.1093/bioinformatics/btab863 Huang D, Su T, Qu L, Wu Y, Gu P He B, Xu X, Zhu C (2016) The complete mitochondrial genome of the Colletes gigas (Hymenoptera, Colletidae, Colletinae). Mitochondrial DNA Part A 27(6): 3878-3879. https://doi.org/10.3109/19401736.2014.987243 Huelsenbeck JP, Ronquist F (2001) MRBAYES: Bayesian inference of phylogenetic trees. Bioinformatics 17(8): 754-755. https://doi.org/10.1093/bioinformatics/17.8.754 Igawa T, Kurabayashi A, Usuki C, Fujii T, Sumida M (2008) Complete mitochondrial ge- nomes of three neobatrachian anurans: a case study of divergence time estimation ZooKeys 1216: 63-82 (2024), DOI: 10.3897/zookeys.1216.131847 79 Jingfang Li et al.: Mitogenomes of two sympatric Rana species using different data and calibration settings. Gene 407(1-2): 116-129. https://doi. org/10.1016/j.gene.2007.10.001 Jiang W, Li J, Xiang H, Sun C, Chang J, Yang J (2023) Comparative analysis and phylo- genetic and evolutionary implications of mitogenomes of Chinese Sinocyclocheilus cavefish (Cypriniformes, Cyprinidae). Zoological Research 44(4): 779-781. https:// doi.org/10.24272/j.issn.2095-8137.2022.439 Lan X, Wang J, Zhang M, Zhou Q, Xiang H, Jiang W (2023) Molecular identification of Acrossocheilus jishouensis (Teleostei, Cyprinidae) and its complete mitochondrial genome. Biochemical Genetics 62(2): 1396-1412. https://doi.org/10.1007/s10528- 023-10501-x Lanfear R, Frandsen PB, Wright AM, Senfeld T, Calcott B (2017) PartitionFinder 2: new methods for selecting partitioned models of evolution for molecular and morpholog- ical phylogenetic analyses. Molecular Biology and Evolution 34(3): 772-773. https:// doi.org/10.1093/molbev/msw260 Lavrov DV, Pett W (2016) Animal mitochondrial DNA as we do not know it: Mt-genome organization and evolution in nonbilaterian lineages. Genome Biology and Evolution 8(9): 2896-2913. https://doi.org/10.1093/gbe/evw195 Levinson G, Gutman GA (1987) Slipped-strand mispairing: a major mechanism for DNA sequence evolution. Molecular Biology and Evolution 4(3): 203-221. https://doi. org/10.1093/oxfordjournals.molbev.a040442 Li J, Chen L, Liu B, Zhao Z, YuL, Yin M, Zhang Y, Gong L (2021) Whole sequence determi- nation and phylogenetic analysis of the mitochondrial genome of the blunt-toothed hairy sea-crab. Journal of Zhejiang Ocean University 40(3): 198-208. Liu C, Hu S (1961) Chinese Tailless Amphibians. Science Publishing House, Beijing, 364 pp. Liu H, Li H, Song F, Gu W, Feng J, Cai W, Shao R (2017) Novel insights into mitochondrial gene rearrangement in thrips (Insecta, Thysanoptera) from the grass thrips, Anapho- thrips obscurus. Scientific Reports 7(1): 4284. https://doi.org/10.1038/s41598-017- 04617-5 Macey JR, Larson A, Ananjeva NB, Fang Z, Papenfuss TJ (1997) Two novel gene orders and the role of light-strand replication in rearrangement of the vertebrate mitoge- nome. Molecular Biology and Evolution 14(1): 91-104. https://doi.org/10.1093/ox- fordjournals.molbev.a025706 Miya M, Nishida M (1999) Organization of the mitogenome of a deep- sea fish Gonosto- ma gracile (Teleostei, Stomiiformes): first example of transfer RNA gene rearrange- ments in bony fishes. Molecular Biology and Evolution 14(1): 416-426. https://doi. org/10.1007/pl00011798 Miya M, Nishida M (2015) The mitogenomic contributions to molecular phylogenetics and evolution of fishes: a 15-year retrospect. Ichthyological Research 62: 29-71. https://doi.org/10.1007/s10228-014-0440-9 Moritz C, Blown WM (1987) Tandem duplications in animal mitochondrial DNAs: varia- tion in incidence and gene content among lizards. Proceedings of the National Acad- emy of Sciences of the United States of America 84(20): 7183-7187. https://doi. org/10.1073/pnas.84.20.7183 Pope CH, Boring AM (1940) A survey of Chinese Amphibia. Peking Natural History Bul- letin 15: 13-86. Rambaut AA (2009) FigTree. Tree Figure Drawing Tool, Version 1.4.2. http://tree.bio. ed.ac.uk/software/Figtree/ [accessed on 6 Jun 2024] Rayko E (1997) Organization, generation and replication of amphimeric genomes: a re- view. Gene 199(1-2): 1-18. https://doi.org/10.1016/S0378-1119(97)00357-0 ZooKeys 1216: 63-82 (2024), DOI: 10.3897/zookeys.1216.131847 80 Jingfang Li et al.: Mitogenomes of two sympatric Rana species Rozas J, Ferrer-Mata A, Sanchez-DelBarrio JC, Guirao-Rico S, Librado P Ramos-On- sins SE, Sanchez-Gracia A (2017) DnaSP 6: DNA sequence polymorphism analysis of large data sets. Molecular Biology and Evolution 34(12): 3299-3302. https://doi. org/10.1093/molbev/msx248 Shao R, Barker SC (2003) The highly rearranged mitogenome of the plague thrips, Thrips imaginis (Insecta, Thysanoptera): convergence of two novel gene boundaries and an extraordinary arrangement of rRNA genes. Molecular Biology and Evolution 20(3): 362-370. https://doi.org/10.1093/molbev/msg045 Shen Y, Jiang J, Yang D (2007) A new species of the genus Rana, Rana hanluica sp. nov. from Hunan Province China (Anura, Ranidae). Acta Zoologica Sinica 53(3): 481-488. Stanton DJ, Daehler LL, Moritz CC, Brown WM (1994) Sequences with the potential to form stem-and-loop structures are associated with coding-region duplications in an- imal mitochondrial DNA. Genetics 137(1): 233-241. https://doi.org/10.1093/genet- leSALS7Z A233 Tan M, Gan H, Lee YP Linton S, Grandjean F (2018) Order within the chaos: insights into phylogenetic relationships within the Anomura (Crustacea, Decapoda) from mito- chondrial sequences and gene order rearrangements. Molecular Phylogenetics and Evolution 127: 320-331. https://doi.org/10.1016/j.ympev.2018.05.015 Tang X, Lyu B, Lu H, Tang J, Meng R, Cai B (2021) The mitochondrial genome of a para- sitic wasp, Chouioiacunea Yang (Hymenoptera, Chalcidoidea, Eulophidae) and phylo- genetic analysis. Mitochondrial DNA Part B 6(3): 872-874. https://doi.org/10.1080/ 23802359.2021.1886008 Van Dijk EL, Auger H, Jaszcezyszyn Y, Thermes C (2014) Ten years of next-generation se- quencing technology. Trends in Genetics 30(9): 418-426. https://doi.org/10.1016/j. tig.2014.07.001 Vinogradov AE, Anatskaya OV (2017) DNA helix: the importance of being AT-rich. Mam- malian Genome 28(9-10): 455-464. https://doi.org/10.1007/s00335-017-9713-8 Wan HL, Yu Z, Qi S, Liu L, Wang Y (2020) A new species of the Rana japonica group (Anura, Ranidae, Rana) from China, with a taxonomic proposal for the R. johnsi group. ZooKeys 942: 141-158. https://doi.org/10.3897/zookeys.942.46928 Wang Y, Shen Y, Feng C, Zhao K, Song Z, Zhang Y, Yang L, He S (2016) Mitogenomic perspectives on the origin of Tibetan loaches and their adaptation to high altitude. Scientific Reports 6: 690. https://doi.org/10.1038/srep29690 Wang Z, Shi X, Guo H, Tang D, Bai Y, Wang Z (2020) Characterization of the complete mitochondrial genome of Uca lacteus and comparison with other Brachyuran crabs. Genomics 112(1): 10-19. https://doi.org/10.1016/j.ygeno.2019.06.004 Wang J, Lan X, Luo Q, Gu Z. Zhou Q, Zhang M, Zhang Y, Jiang W (2022) Characteriza- tion, comparison of two new mitogenomes of crocodile newts Tylototriton (Cauda- ta, Salamandridae), and phylogenetic implications. Genes 13(10): 1878. https://doi. org/10.3390/genes13101878 Wu Z (1981) Karyotype of Chinese wood frogs from Beijing. Journal of Genetics 8: 138-144. Xia X (2018) DAMBE7: New and improved tools for data analysis in molecular biolo- gy and evolution. Molecular Biology and Evolution 35(6): 1550-1552. https://doi. org/10.1093/molbev/msy073 Xia X, Li Y, Yang D, Pi Y (2022) Habitat characteristics during the breeding season of the cold-dew wood frog (Rana pipiens) and the main factors affecting its habitat selec- tion. Journal of Ecology 41(9): 1740-1745. Xiang H, Zhou Q, Li W, Shu J, Gu Z, Jiang W (2024) Insights into phylogenetic positions and distribution patterns: Complete mitogenomes of two sympatric Asian horned ZooKeys 1216: 63-82 (2024), DOI: 10.3897/zookeys.1216.131847 8] Jingfang Li et al.: Mitogenomes of two sympatric Rana species toads in Boulenophrys (Anura: Megophryidae). Ecology and Evolution 14(7): e11687. https://doi.org/10.1002/ece3.11687 Yan F, Jiang K, Chen H, Fang P Jin J, Yi Li, Wang S, Murphy RW, Che J, Zhang Y (2011) Matrilineal history of the Rana longicrus species group (Rana, Ranidae, Anura) and the description of a new species from Hunan, southern China. Asian Herpetological Research 2(2): 61-71. https://doi.org/10.3724/SP.J.1245.2011.00061 Yang H, Li T, Dang K, Bu W (2018) Compositional and mutational rate heterogeneity in mitochondrial genomes and its effect on the phylogenetic inferences of Cimicomor- pha (Hemiptera, Heteroptera). BMC Genomics 19(1): 264. https://doi.org/10.1186/ $12864-018-4650-9 Ye F, Samuels DC, Clark T, Guo Y (2014) High-throughput sequencing in mitochondrial DNA research. Mitochondrion 17: 157-163. https://doi.org/10.1016/j.mito.2014.05.004 Ye N, Wang X, Li J, Bi C, Xu Y, Wu D, Ye Q (2017) Assembly and comparative analysis of complete mitochondrial genome sequence of an economic plant Salix suchowensis. PeerJ 5: e3148. https://doi.org/10.1186/s12864-021-07490-9 Ye F, Li H, Xie Q (2021) Mitogenomes from two specialized subfamilies of Reduviidae (Insecta, Hemiptera) reveal novel gene rearrangements of true bugs. Genes 12(8): 1134. https://doi.org/10.3390/genes1 2081134 Zhang Y, Gong L, Lu X, Jiang L, Liu B, Liu L, Lu Z, Li RP Zhang X (2020) Gene rearrange- ments in the mitochondrial genome of Chiromantes eulimene (Brachyura, Sesarmi- dae) and phylogenetic implications for Brachyura. International Journal of Biological Macromolecules 162: 704-714. https://doi.org/10.1016/j.ijbiomac.2020.06.196 Zhang M, Zhou Q, Xiang H, Wang J, Lan X, Luo Q, Jiang W (2023) Complete mitochon- drial genome of Rectoris luxiensis (Teleostei, Cyprinidae): characterization and phylo- genetic implications. Biodiversity Data Journal 11: e96066. https://doi.org/10.3897/ BDJ.11.e96066 Zhi M , Li W (2016) Comparison of skin histostructures in three species of genus Rana. Chinese Journal of Zoology 51(5): 844-852. Zhou K, Li H, Han D, Bauer AM, Feng J (2006) The complete mitochondrial genome of Gekko gecko (Reptilia: Gekkonidae) and support for the monophyly of Sauria includ- ing Amphisbaenia. Molecular Phylogenetics and Evolution 40(3): 887-892. https:// doi.org/10.1016/j.ympev.2006.04.003 Zhou Y, Wang §S, Zhu H, Li PR Yang B, Ma J (2017) Phylogeny and biogeography of South Chinese brown frogs (Ranidae, Anura). PLOS ONE 12(4): e175113. https://doi. org/10.1371/journal.pone.0175113 Zhou Q, Xiang H, Zhang M, Liu Y, Gu Z, Lan X, Wang J, Jiang W (2023) Two complete mitochondrial genomes of Leptobrachium (Anura, Megophryidae, Leptobrachiinae): characteristics, population divergences, and phylogenetic implications. Genes 14(3): 768. https://doi.org/10.3390/genes14030768 Zhu S, Zhang X, Zhang S, Zhang R, Liu S, Xiong R (2018) Molecular identification of a Rana hanluica population distributed in Youyang County, Chongqing City in China. Open Journal of Nature Science 6(1): 1-8. https://doi.org/10.12677/ojns.2018.61001 Zyla J, Marczyk M, Domaszewska T, Kaufmann SHE, Polanska J, Weiner J (2019) Gene set enrichment for reproducible science: comparison of CERNO and eight other algo- rithms. Bioinformatics 35(24): 5146-5154. https://doi.org/10.1093/bioinformatics/ btz447 ZooKeys 1216: 63-82 (2024), DOI: 10.3897/zookeys.1216.131847 99