Phylogeny and evolution of Lasiopodomys in subfamily Arvivolinae based on mitochondrial genomics

The species of Lasiopodomys Lataste 1887 with their related genera remains undetermined owing to inconsistent morphological characteristics and molecular phylogeny. To investigate the phylogenetic relationship and speciation among species of the genus Lasiopodomys, we sequenced and annotated the whole mitochondrial genomes of three individual species, namely Lasiopodomys brandtii Radde 1861, L. mandarinus Milne-Edwards 1871, and Neodon (Lasiopodomys) fuscus Büchner 1889. The nucleotide sequences of the circular mitogenomes were identical for each individual species of L. brandtii, L. mandarinus, and N. fuscus. Each species contained 13 protein-coding genes (PCGs), 22 transfer RNAs, and 2 ribosomal RNAs, with mitochondrial genome lengths of 16,557 bp, 16,562 bp, and 16,324 bp, respectively. The mitogenomes and PCGs showed positive AT skew and negative GC skew. Mitogenomic phylogenetic analyses suggested that L. brandtii, L. mandarinus, and L. gregalis Pallas 1779 belong to the genus Lasiopodomys, whereas N. fuscus belongs to the genus Neodon grouped with N. irene. Lasiopodomys showed the closest relationship with Microtus fortis Büchner 1889 and M. kikuchii Kuroda 1920, which are considered as the paraphyletic species of genera Microtus. TMRCA and niche model analysis revealed that Lasiopodomys may have first appeared during the early Pleistocene epoch. Further, L. gregalis separated from others over 1.53 million years ago (Ma) and then diverged into L. brandtii and L. mandarinus 0.76 Ma. The relative contribution of climatic fluctuations to speciation and selection in this group requires further research.

In the present study, we sequenced the whole mitochondrial genomes of L. mandarinus, L. brandtii, and N. fuscus, which are species with three repeat individuals, using highthroughput sequencing technology and used the complete mitochondrial genomes of related species from the National Center for Biotechnology Information database to clarify the generic taxonomy of Lasiopodomys and evolutionary history of adaptation on aboveground and subsurface life. The findings of this research provide evolutionary information regarding the hypoxia adaptation of Lasiopodomys.
The phylogenetic relationships of the two different matrices as well as the whole mitochondrial genomes and PCG sequence matrices were constructed using the maximum likelihood (ML) approach in IQ-TREE v1.6.12 (Nguyen et al., 2015) and Bayesian analysis (BI) in the BEAST v1.8.4 program (Drummond & Rambaut, 2007). We conducted analysis using 5000 ultrafast bootstrap replicates and the best-fit model in the IQ-TREE software.
To determine the maximum clade credibility trees of two different matrices, BEAST analyses were performed using the GTR+G+I substitution models identified above and the uncorrelated relaxed clocks for clock type (Drummond et al., 2006), Yule process for tree prior (Gernhard, 2008), and other default parameters. Each Markov chain Monte Carlo of 20,000,000 generations was sampled in every 10,000 generations. The effective sample sizes were estimated using Tracer v1.7 for all parameters more than 200 (Rambaut et al., 2018). Maximum clade credibility trees were constructed using TreeAnnotator v1.8.4 with a burn-in of the first 20% of the sampled trees (Drummond & Rambaut, 2007). Positive selection in all 13 PCGs was determined using branch models and branch-site models via phylogenetic analysis using ML (PAML4.7) programs (Yang, 2007). Branch models were used with the one-ratio model, i.e., all the species had the same ω ratio, and the ω = 1 model, with all species in natural selection. Based on the phylogenetic tree, we estimated the ω values of each PCG. The branch-site models used all Lasiopodomys species as the foreground branches, and the likelihood ratio test (LRT) was conducted to assess the statistical significance of positive selection.
The molecular divergence time was estimated using the Yule and birth-death processes for trees before implementing phylogeny construction using BEAST v1.8.4 (Gernhard, 2008;Heath, Huelsenbeck & Stadler, 2014). Marginal likelihood estimation for path sampling and stepping-stone sampling (Xie et al., 2011) using 5,000,000 in chain lengths of 500 path steps was used to sample the likelihood of every 5,000 chains (Baele et al., 2012;Baele et al., 2013). We applied three constraints to calibrate the tree at three prior nodes: (1) the divergence time of the Taiwan vole, Microtus kikuchii Kuroda 1920, and the reed vole Microtus fortis, of which the split between the subgenus Alexandromys Ognev 1914 and Pallasiimus Schrank 1798 was estimated via molecular clock analysis at ∼1.19 ± 0.19 Ma (Bannikova et al., 2010;Gao et al., 2017), (2) the earliest known fossil of Eothenomys Allen 1924 at 2.0 Ma (Liu et al., 2012a;Kohli et al., 2014), and (3) the oldest fossil of Arvicola, which was estimated at 3.0-3.5 Ma (Abramson et al., 2009a;Chen et al., 2012); we used the mean value of 3.25 Ma.

Ecological niche modeling
The maximum entropy (Maxent) method was used to predict the current potential geographic distributions of L. mandarinus, L. brandtii, L. gregalis, and N. fuscus as well as their suitable distributions during the mid-Holocene, 6,000 years ago (kya), Last Glacial Maximum (LGM; 22 kya), and Last Interglacial (LIG; 120-140 kya) epochs (Phillips, Anderson & Schapire, 2006;Elith et al., 2011). Presence records were obtained for all four species according to the GBIF database and published papers (Appendix S2). Climatic variables with 19 bioclimatic layers were obtained from the database WorldClim version 1.4 at a resolution of 2.5 arc-minute grid format (Hijmans et al., 2005). The potential distributions of the species during the LGM and Holocene periods were predicted using both MIROC-ESM and CCSM4 models (Watanabe et al., 2011;Shields et al., 2012). Strongly correlated bioclimatic layers (r > 0.9) as determined using Pearson's correlation analysis in R 3.6.2 (Appendix S3) (R Development Core Team, 2013) were excluded. Moreover, Maxent was independently performed among these species using area under the receiver operating characteristic curve (AUC) prediction model evaluation (DeLong, DeLong & Clarke-Pearson, 1988;Fawcett, 2006).

RESULTS
The whole mitochondrial genome length of L. mandarinus was 16,562 bp, with the same sequences among repeated individuals. The mitochondrial genome length of L. brandtii was only 5 bp shorter than that of L. mandarinus, whereas that of N. fuscus was 220 bp shorter than that of L. mandarinus (Fig. 1). On the other hand, L. mandarinus was found to be 234 bp longer than the former sequenced mitogenomes (GenBank no. KF819832 & JX014233). All sequences of the three species were longer than those of L. gregalis, a species previously in the genus Microtus, with sequence lengths of 16,292 bp (GenBank no. MN199169) and 16,294 bp (GenBank no. MN199170). All the three mitogenomes were assembled into a typical circular map with 13 PCGs, 22 tRNAs, 2 rRNAs (rrn12 and rrn16), and a D-loop region (Fig. 1, Table 1). Five types of start codons-ATA, ATC, ATG, ATT, and GTG-were identified among the PCGs, whereas three types of stop codons were identified for these species. The nucleotide composition of L. brandtii, L. mandarinus, and N. fuscus was biased for A+T by 59.5%, 59.5%, and 58.4%, respectively. All these mitogenomes showed a positive AT skew of 0.08 for L. brandtii, 0.09 for L. mandarinus, and 0.09 for N. fuscus. However, these species showed a negative GC skew ranging from −0.30 for L. brandtii to −0.34 for L. mandarinus ( Fig. 1, Table 2). L. gregalis showed higher AT skew (0.10) and GC skew (−0.30) compared with the other three species. Among the 13 PCGs in these 4 species, nucleotide composition ranged from −0.69 in ATP8 to −0.16 in ND4L for L. mandarinus, with a GC skew ranging from −0.14 in ND4L for L. brandtii to 0.33 in ND6 for L. mandarinus. Similarly, all 13 PCGs exhibited a negative GC skew; however, COX1, ND4L in all species, COX3 in L. brandtii and L. mandarinus, and ND3 in N. fuscus showed a negative AT skew and ND3 in L. brandtii and L. mandarinus had an AT skew of 0 ( Table 2).
The nucleotide diversity among the published Arvicolinae mitogenome sequences and our study species was 0.1429 ± 0.0001, whereas the nucleotide diversity of the mitogenomes of Lasiopodomys was 0.0836 ± 0.0155 (Fig. 2). The total nucleotide diversity in all 13 PCGs of Arvicolinae and the genus Lasiopodomys was 0.1603 ± 0.0027 and 0.0953 ± 0.0180, respectively (Fig. 2). In Arvicolinae, nucleotide diversity ranged from 0.1378 ± 0.0049 in Cytb to 0.1977 ± 0.0077 in ND3, whereas for Lasiopodomys, it ranged from 0.0829 ± 0.0157 in COX3 to 0.1256 ± 0.021 in ND4L. The results of the ML and Bayesian approaches were applied to the datasets of the whole mitogenomes, and the 13 PCG matrices inferred the same topology of the phylogenetic tree structure (Fig. 3). Our results supported that Lasiopodomys, Microtus, and Neodon have close relationships with the basal group of Proedromys Thomas 1911. Furthermore, the phylogenetic tree suggested that L. brandtii, L. mandarinus, and L. gregalis formed the genus of Lasiopodomys, whereas N. fuscus showed a close relationship with N. irene, belonging to the genus Neodon. Microtus was subdivided into two groups: one containing M. fortis and M. kikuchii, which were strongly supported as the sister group to Lasiopodomys, and the other was the basal group of the above species.
In the branch models, the one-ratio model was determined as superior to the ω = 1 model (df = 1, p < 0.01), suggesting that all the PCGs in the mitogenomes of Lasiopodomys undergo purifying selection (Table 3). In the branch-site model, only the ATP6 gene was present in some positive selection sites (60I 0.987, p < 0.01) in Lasiopodomys (Table 3). Moreover, positive selection sites were predicted in Cox1, Cox3, Cytb, ND2, ND3, and ND5. However, the LRTs were not significant.   The species divergence time among the Lasiopodomys species and related genera was calculated using the uncorrelated relaxed molecular clock model, which was calibrated with three prior divergence times of Arvicolinae (Fig. 3). The results suggested that the origin of Lasiopodomys was no earlier than the early Pleistocene epoch (∼0  The high AUC values determined via ecological niche modeling (ENM) indicated the good performance of the model predictions of this study (Appendix S4). During the periods from the LIG to present, all species of Lasiopodomys showed no evident loss of a suitable habitat. A western expansion of L. brandtii has been predicted in Northeast China, Inner Mongolia, and South Siberia, whereas a weak fragment was predicted for L. gregalis among the Eurasia regions (Fig. 4). Moreover, suitable areas were predicted in highly suitable habitat regions during the LGM in these species. More northern suitable areas were predicted during the LIG, and a northern expansion was predicted during the transition from the Holocene period to the present (Fig. 4). In addition, highly suitable habitats were observed for N. fuscus in the Hengduan Mountains during all periods, whereas more eastern distributions were predicted during the LGM (Fig. 4).

Structural features of the whole mitochondrial genome of Lasiopodomys
Among the nine complete mitochondrial sequences, all the species showed same sequences in the three repeated individuals, thereby supporting the accuracy and low intraspecific variation of our studies (Brown & Simpson, 1981). Although N. fuscus showed similar characteristics to previously sequenced mitogenomes (GenBank no. MG833880), L. mandarinus exhibited a longer sequence than that previously reported (Cong et al., 2016;Li, Lu & Wang, 2016a;Li et al., 2016b;Li et al., 2019). This difference may be due to nucleotide errors, particularly in tandem repeats, caused by different sequencing technologies: Sanger sequencing versus high-throughput sequencing (Pfeiffer et al., 2018). All these differences occurred in the intergenic region, with little impact on subsequent analysis. Therefore, we reserved both types of sequence data in the subsequent analysis.
All the PCGs of these species, similar to the other Arvicolinae mitogenomes, had an incomplete stop codon that was automatically filled during the transcription process in the mitogenomes of animals, with no effect on translation (Ojala, Montoya & Attardi, 1981). Similar to previous studies, the nucleotide diversity of all the PCGs in both Lasiopodomys and Arvicolinae typically showed the highest divergence in the NADH dehydrogenase complex and the lowest divergence in the cytochrome c oxidase subunit complex and cytochrome B gene (Huang et al., 2019;Ramos et al., 2018). The nucleotide sequence diversity of the NADH dehydrogenase gene groups may be affected by variations in the historical environment (Ramos et al., 2018;Mueller, 2006). Similar to previously published mitogenomes, the AT skew of Lasiopodomys and N. fuscus was consistent with that of vertebrates (Zhang, Cheng & Ge, 2019;Martin, 1995), further indicating evolutionary pressure related to the mechanism of DNA replication (Charneski et al., 2011;Dai & Holland, 2019).

Phylogenetic relationships of Lasiopodomys
Our molecular phylogenetic analysis results were highly consistent those of previous studies. In our study, the subfamily Arvicolinae was supported as a monophyletic group based on the molecular data of Cytb, COX1, GHR, CLOCK, and BMAL1 (Abramson et al., 2009b;Buzan et al., 2008;Liu et al., 2017;Martin et al., 2000;Sun et al., 2018). Our results suggest that N. (Lasiopodomys) fuscus within the genus Neodon forms a sister relationship with N. irene, consistent with the results reported by Chen et al. (2012) andLi et al. (2019). The stable clustering of L. brandtii, L. mandarinus, and L. gregalis into one group confirms the systematic positions of Lasiopodomys. This topology was consistent with that of other phylogenetic studies based on nuclear genes (Sun et al., 2018), mitochondrial DNA (Abramson et al., 2009a;Liu et al., 2012b;Martínková & Moravec, 2012;Petrova et al., 2016), and whole genomes (Li, Lu & Wang, 2016a;Li et al., 2019;Tian et al., 2020). However, it contradicts with the systematic position based on the morphological characteristics of these species (Allen, 1940;Corbet, 1978;Wilson & Reeder , 2005). Further, L. brandtii and L. mandarinus have consistently presented as a sister group in molecular phylogenetic studies, with seldom distinguished morphological characteristics but different aboveground and underground habitats, suggesting a mechanism of environmental adaptation during rapid speciation (Alexeeva, Erbajeva & Khenzykhenova, 2015;Dong et al., 2018;Li et al., 2017). Other species of Microtus and Neodon were not found in the monophyletic group (Liu et al., 2012a); M. kikuchii and M. fortis were grouped as sister lineages within the Lasiopodomys clades and were considered belonging to the subgenus Alexandromys based on phylogenetic research (Mezhzherin, Zykov &Morozov-Leonov, 1993), allozymes, andCytb (Bannikova et al., 2010). All these genera form a ''Microtus s. l.,'' which could be the ''core Arvicolinae' ' (Baca et al., 2019).

Evolution and demographic history of Lasiopodomys
When inferring the divergence time of Lasiopodomys and related genera, both the Yule process and birth-death process speciation models were required with multiple fossil calibration nodes employed in phylogenetic analysis to develop more robust estimates (Drummond & Rambaut, 2007;Humphreys et al., 2016). Based on complete genomes and PCG phylogenetic trees, both models presented similar estimates of a relatively recent origin and divergence time for Microtus s. l. during the early Pleistocene epoch. The oldest reported fossil of Microtus s. l. was during the early Pleistocene epoch (Chaline et al., 1999). An arid and cold environment raised species dispersal and speciation in response to Pleistocene climatic fluctuations (Vasconcellos et al., 2019). Our study supported the first appearance of Lasiopodomys in the late early Pleistocene epoch from the Transbaikal area (Alexeeva, Erbajeva & Khenzykhenova, 2015;Li et al., 2017) at ∼1.52-2.09 Ma (Petrova et al., 2016) but later than that estimated by chromosomes at 3 Ma (Gladkikh et al., 2016). At ∼1.28-1.81 Ma, the morphological characters of L. gregalis proposed the earliest clades of modern Lasiopodomys, as indicated by molecular data and fossils (Abramson et al., 2009a;Chaline et al., 1999;Petrova et al., 2016). Thereafter, the clades separated into L. brandtii and L. mandarinus at ∼0.58-0.98 Ma in our study, which is similar to inferences from Cytb and D-loop sequences (Li et al., 2017;Petrova et al., 2015) but less similar to the inferences from molecular cytogenetic analyses at ∼1.8 Ma (Gladkikh et al., 2016).
ENM indicated a considerably wider distribution area of Lasiopodomys in the past than in the present, which conforms to the fossils from the Pleistocene period (Alexeeva, Erbajeva & Khenzykhenova, 2015). During the early Pleistocene period, continuous cooling formed an arid climate in the high latitudes of the Northern Hemisphere (Guo et al., 2008). Climatic changes seldom shifted the suitable habitat of Lasiopodomys during the LIG and LGM periods. It is possible to infer that migration events occurred during the extremely cold and dry conditions, with a trend of continuous distribution farther to the northeast during the Pleistocene period until the Holocene period (Alexeeva, Erbajeva & Khenzykhenova, 2015;Prost et al., 2013). The appearance of N. fuscus, which is adapted to plateau climates, was later than the Qinghai-Tibet Plateau uplift (Wang et al., 2008), with no significant distributed shifts. All ancient species of Lasiopodomys may have been distributed as per their current distribution areas with a radiation evolution (Abramson et al., 2009b;Bannikova et al., 2010) before the interglacial and glacial periods based on ENM and fossil reports (Alexeeva, Erbajeva & Khenzykhenova, 2015;Petrova et al., 2015). Considering the lower sensitivity to climatic changes and adaptation to habitat areas, the Lasiopodomys species could colonize in north regions; moreover, the evolution of characteristics, such as teeth and densely furred plantar surfaces, further enabled their survival in cooler, drier conditions.
Despite precipitation and temperature fluctuations, a decline in atmospheric O 2 also occurred during the past 0.8 Ma (Stolper et al., 2016). Environmental stress caused a major driving on evolutionary process (Parsons, 2005). In the species of rodents, limited oxygen availability resulted in evolutionary adaptation and appearance of various strategies (Pamenter et al., 2020), such as different colonial habitats of life-subterranean (L. mandarinus) and plateau (L. fuscus); these strategies formed unique physiological and molecular adaptations to hypoxia (Jiang et al., 2020;Dong, Wang & Jiang, 2020). Our study supports a history of rapid population expansion under positive selection via mitogenome sequences such as the ATP6 gene, which uses oxygen to create adenosine triphosphate. However, further research using integrated phylogeographic analyses of the genus Lasiopodomys (Li et al., 2017;Petrova et al., 2015) is warranted to determine the adaptation of L. brandtii and L. mandarinus to factors including precipitation, temperature, and chronic hypoxia.