The phylogeny of the Anderson's White‐bellied Rat (Niviventer andersoni) based on complete mitochondrial genomes

Abstract The phylogenetic structure of the genus Niviventer has been studied based on several individual mitochondrial and nuclear genes, but the results seem to be inconsistent. In order to clarify the phylogeny of Niviventer, we sequenced the complete mitochondrial genome of white‐bellied rat (Niviventer andersoni of the family Muridae) by next‐generation sequencing. The 16,291 bp mitochondrial genome consists of 22 transfer RNA genes, 13 protein‐coding genes (PCGs), two ribosomal RNA genes, and one noncoding control region (D‐Loop). Phylogenetic analyses of the nucleotide sequences of all 13 PCGs, PCGs minus ND6, and the entire mitogenome sequence except for the D‐loop revealed well‐resolved topologies supporting that N. andersoni was clustered with N. excelsior forming a sister division with N. confucianus, which statistically rejected the hypothesis based on the tree of cytochrome b (cytb) gene that N. confucianus is sister to N. fulvescens. Our research provides the first annotated complete mitochondrial genome of N. andersoni, extending the understanding about taxonomy and mitogenomic evolution of the genus Niviventer.

) (Ge et al., 2017(Ge et al., , 2021. They live in various kinds of forests in both lowlands and mountains (Li & Yang, 2009). Fossil records showed that this species extended to the low-altitude regions of Southeast China during the late Quaternary in Chongqing and Guizhou, suggesting they migrated southward when the climate turned colder (Bahain, 2007;Bekken et al., 2004).
The phylogenetic position of N. andersoni has not been fully determined, due to controversial phylogenetic topologies within the genus Niviventer. An early study proposed that the Niviventer could be divided into two primary groups: the N. andersoni-Division and the N. niviventer-Division (Musser, 1981). Phylogenetic trees based on mitochondrial cytochrome b (cytb) gene showed that N. andersoni and N. excelsior were clustered together and comprised the N.
andersoni-Division (Jing et al., 2007). Meanwhile, N. confucianus, together with N. fulvescens and N. cremoriventer, formed N. niviventer-Division, a sister division to N. andersoni-Division ( Figure 1) (Ge et al., 2021;He & Jiang, 2015;Lu et al., 2015;Zhang et al., 2016). However, this classification might be inconclusive due to the limited number of gene sequences used, since each gene evolves under different evolutionary pressures and time scales (Choi & Kim, 2017). Compared to individual mitochondrial gene sequences, complete mitochondrial genome sequences can provide higher resolution and sensitivity for better revealing the evolutionary relationships among closely related species (Ladoukakis & Zouros, 2017;Wei et al., 2017).
Up to now, the mitogenomes of the four species (N. andersoni, N. confucianus, N. fulvescens, and N. niviventer) within the genus Niviventer have not been entirely sequenced or verified. Since complete mitochondrial genomes have been used for taxonomic and phylogenetic analyses of diverse animal groups (Kim et al., 2020;Ladoukakis & Zouros, 2017;Lavrov & Pett, 2016;Wei et al., 2017), the lack of genetic data has limited our understanding of the phylogeny of N. andersoni. In the present study, we sequenced the complete mitochondrial genome of N. andersoni and described typical features of the N. andersoni mitochondrial genome revealing its phylogenetic relationships with other white-bellied rat species. Our findings highlight the importance of complete mitogenome information in phylogenetic analyses for rodent species.

| Sample collection and genomic DNA extraction
Individuals of N. andersoni were collected from Lufeng County, Yunnan province, China (H = 1875.43 m), in August 2018. These individuals were sacrificed and dissected for organ collection. The heart, liver, spleen, lung, kidney, and muscle were kept in the cryopreservation tubes directly. All the samples were immediately put in liquid nitrogen for short storage, then transported to the laboratory in dry ice and stored at −80°C. DNA was extracted from the muscle using mitochondrial extraction kit (Solarbio) and stored at −80°C.

| Mitogenome sequencing, assembly, and annotation
The mitochondrial DNA was subjected to random PCR (rPCR) as previously described (Li et al., 2010). Briefly, the extracted mitochondrial DNA was first synthesized by random primer (5'-GCC GGA GCT CTG CAG AAT TCNNNNNN-3'). Then 2 μl of the synthetic product was used to perform the PCR in a 50μl reaction mixture volume containing 10 μl PCR buffer, 1 mM MgSO4, 0.2 mM each dNTP, 40 pmol universal primer (5'-GCC GGA GCT CTG CAG AAT TC-3'), F I G U R E 1 Niviventer andersoni. This photograph is adapted from Ge et al. (2021), Molecular Biology and Evolution, 2021 (Wilson & Reeder, 2007) under proper copyright permission. The schematic diagram of the evolutionary tree embedded in the bottom of the picture shows the N. niviventer-Division and N. andersoni-Division and 1 U KOD-Plus DNA polymerase (Toyobo, Japan). The reaction was conducted for 40 cycles of 94°C for 30 s, 54°C for 30 s, and 68°C for 2 min followed by incubation for 10 min at 68°C. The products were analyzed by agarose gel electrophoresis.
The purified rPCR products were used to construct the sequencing library and sequenced on HiSeq-PE150 instrument (TIANGEN, Beijing, China). The raw reads were trimmed and filtered using Trimmomatic (Version 0.39) (Bolger et al., 2014). The cleaned reads were aligned to NCBI non-redundant protein sequence database using BLASTx by DIAMOND (Buchfink et al., 2015). Mitochondrial reads were selected and de novo assembled into a complete mitochondrial genome using Geneious software package (Version 2019.1.1) (Kearse et al., 2012).
Protein-coding genes (PCGs) were annotated using the NCBI ORF Finder (https://www.ncbi.nlm.nih.gov/orffi nder/) and BLASTx with the vertebrate mitochondrial genetic code. The tRNA genes were identified using the tRNAscan-SE Search Server under the default search mode, using the vertebrate mitochondrial genetic code source (Chan & Lowe, 2019). Composition skew analysis was calculated according to the formulas: AT skew = (A-T)/(A + T) and GC skew = (G-C)/(G + C) (Perna & Kocher, 1995). Relative synonymous codon usage (RSCU) values were calculated using CodonW 1.4.2 . The circular mitochondrial genome map of N. andersoni was drawn using OGDRAW   (Table 1).

| Phylogenetic analysis
Meanwhile, cytb, cox1, and ND2, the three genes commonly for phylogenetic analyses, were used to construct phylogenetic control trees.
To further investigate the phylogenetic relationships of N. andersoni, the phylogenetic relationships were reconstructed based on the complete mitochondrial genome. The D-loop region was excluded because of the rapid mutation rate in this region. The nucleotide sequences were aligned using ClustalX with default settings before concatenation by DAMBE (Version 7.2) (Larkin et al., 2007;Xia, 2017). Models of evolution were evaluated using corrected Aikake Information Criteria (AICc) in jModelTest 2.1.10 to determine the best nucleotide substitution model (Darriba et al., 2012). Maximum likelihood (ML) analysis of the 13 PCGs in 13 species of rodent was also performed using MEGA X (Kumar et al., 2018). The support values of the ML tree were evaluated via a bootstrap test with 1000 iterations.

| Estimation of divergence date
The 13 PCG sequences were aligned using Muscle program by codon method in MEGA. We calculated the differentiation time of species using 13 PCGs as different partitions in BEAST v1.10.4 program . The substitution models and clock models (uncorrelated lognormal relaxed clock) were unlinked among partitions in BEAST. The most appropriate substitution and variant sites model of each PCG through ModelFinder program was estimated according to the Bayesian Information Criterion (BIC) method (Kalyaanamoorthy et al., 2017), the result is shown in Table S1. Referring to previous reports of divergence dates in Rattus genus, two priors of the most recent common ancestor (tMRCA) were used to calibrate the molecular clock (Camacho-Sanchez & Leonard, 2020;Robins et al., 2008). The

| Genome organization
From the raw reads, a total of 1,578,672 high-quality reads were obtained and used to assemble the N. andersoni mitochondrial genome.
As a result, the complete mitochondrial genome sequence of N. andersoni was deposited into NCBI with GenBank accession number MW030174. The mitogenome of N. andersoni was a circular DNA molecule with 16,291 bp in length. As shown in Figure 2, the mitogenome organization of N. andersoni was similar to those of other rodents (Boore, 1999). Thirty-seven typical mitochondrial genes were identified, including 13 PCGs, 22 tRNA genes, and 2 rRNA genes (Table S2).
Most of the genes were encoded by the heavy strand (H-strand), while ND6 and 8 tRNAs were encoded by the light strand (L-strand).
The total base composition of N. andersoni mitochondrial genome was estimated to be 33.7% for A, 25.8% for C, 12.1% for G, and 30.0% for T, making AT and GC percentage as 61.6% and 38.4%, respectively, indicating that the mitochondrial genome biased toward AT (Table 2). Such base composition bias has been reported to play a vital role in the replication and transcription of mitochondrial genome (Wei et al., 2010). It also showed a negative GC skew value (−0.347), indicating that C was more common than G whereas the AT skewness was positive (0.092) suggesting that A occurred more frequently than T in the N. andersoni mitochondrial genome (Table 2).

| Protein-coding genes
The total length of the 13 PCGs was 11,420 bp, composing 70.

| Ribosomal RNA and Transfer RNA genes
The mitogenome of N. andersoni contained the typical 22 tRNA genes throughout the genome and appeared to be highly A+T biased, ranging in length from 59 bp to 75 bp. Among these tRNA genes, 8 were transcribed on the L-strand and the remaining 14 were transcribed on the H-strand (Table S2). All the tRNA genes exhibited a typical cloverleaf structure, except trns1 which lacked a dihydroxyuridine (DHU) arm that was simplified to a ring shape. Loss of the DHU arm is common in the mitogenomes of many mammal animals (Wolstenholme, 1992).
The two rRNA genes (lrRNA, srRNA) encoding the small and large ribosomal subunits were located between tRNA Phe and tRNA Leu on the L-strand of N. andersoni. The lrRNA and srRNA genes were 1567 and 957 bp in length, respectively. The A+T content of rRNA was 63.43%, and its AT skew (0.204) and GC skew (−0.099) showed that more As and Cs were present in the rRNA than As and Gs (Table 2).

| Phylogenetic analysis
Based on 13 PCGs of 13 species, we obtained a phylogenetic tree by ML method with 1000 replications in which Mus musculus was set as the outgroup (Figure 4a). Previous research has suggested that ND6 gene should be excluded during phylogenetic analysis due to its high heterogeneity and consistently poor phylogenetic performance (Miya & Nishida, 2000). Thus, we constructed another phylogenetic tree based on PCGs excluding ND6 (Figure 4b). As a result, the two phylogenetic analyses were similar. When compared with other white-bellied rat species, N. andersoni was phylogenetically closer to N. excelsior and clustered within genus Niviventer. For the three independent trees generated for cytb, cox1, and ND2, we obtained three different topological structures. The ND2 tree showed similar topologies to the 13 PCGs combined, both showing that N. andersoni and N. confucianus were sister species. The cox1 tree showed a cluster of N. confucianus, N. cremoriventer, and N. fulvescens. For cytb, N. confucianus was clustered with M. musculus ( Figure S1).
The ML tree constructed based on the complete mitochondrial genome (except D-loop) showed the same topologies as those PCG trees ( Figure 5). In addition, we estimated that the divergence date

| DISCUSS ION
In this study, we obtained the first annotated complete mitochon- The phylogenetic topologies within the genus Niviventer have been studied for decades. This genus was initially divided into two groups based on cytb gene that were the N. andersoni-Division  . andersoni and N. excelsior, and the N. niviventer-Division including N. confucianus, N. fulvescens, and N. cremoriventer (Jing et al., 2007). However, as more Niviventer species were found and phylogeny on several other mitochondrial and nuclear genes was     (Ge et al., 2021;He & Jiang, 2015;Lu et al., 2015;Zhang et al., 2016).
Although the latter classification has been accepted by most researchers recently, the phylogenetic topologies were still based on a limited number of genes only, which might lead to bias. Therefore, more genetic sequence, especially the nuclear genome, is required to further clarify it.
Our results based on the complete mitochondrial genome partially supported the latter four-division taxa. The results suggest that N. andersoni and N. excelsior clustered together, then with N.
confucianus, and these three formed a division, and N. fulvescens and N. cremoriventer formed another division. Since each gene evolves under different evolutionary pressure and time scale, it has been known that one gene tree for a population may differ from other gene trees for the same population depending on the subjective selection of the genes (Choi & Kim, 2017). Our results also suggested that the inconsistency among evolutionary trees of mitochondrial single gene could be verified by mitochondrial whole genome evolutionary analysis.
Nevertheless, the phylogenetic topologies of the genus Niviventer revealed by maternally inherited mitochondrial genome alone may not agree with that obtained with nuclear genome, due to incomplete lineage sorting, mitochondrial introgression, and recent hybridization (Barbosa et al., 2018;Rubinoff & Holland, 2005 Therefore, more nuclear genes or even complete nuclear genome should be included in the analysis to provide more comprehensive view of the phylogenetic topologies in the genus Niviventer.

| CON CLUS ION
We have sequenced and annotated the complete mitochondrial

ACK N OWLED G M ENTS
We thank Dr. Deyan Ge from Key Laboratory of Zoological Systematics and Evolution, Institute of Zoology, Chinese Academy

F I G U R E 6
The tree with divergence time based on 13 PCGs. The value near the node indicates the age of the node, and the node shape indicates the 95% HPD range. The label on the branch represents the Bayesian posterior probability of Sciences for providing us with the photograph of Niviventer ander-