Although taxonomical and molecular systematics have led to some progress in the relationship between the genus Lasiopodomys and its related genera, numerous uncertainties remain unelucidated. The species belonging to this genus was first described by Lataste in 1887 as part of the Arvivolinae Gray 1821 (Cricetidae Fischer 1817) subfamily, which includes the genera Phaiomys Blyth 1863, Microtus Schrank 1798, and Neodon Horsfield 1841 (Allen, 1940; Corbet, 1978; Gromov & Polyakov, 1978; Liu et al., 2013; Wang, 2003). The genus Lasiopodomys includes three species from different colonial habitats of life—subterranean (L. mandarinus Milne-Edwards 1871), aboveground (L. brandtii Radde 1861), and plateau (L. fuscus Büchner 1889) by Wilson & Reeder (2005)—with relatively short tail and densely furred plantar surfaces. However, their generic taxonomy is not universally accepted, specifically in relation to Phaiomys, Microtus, and Neodon. Molecular data have revealed that the narrow-headed vole Microtus gregalis Pallas 1779 (formerly included in subgenus Stenocranius Katschenko 1901) is closely related to the species belonging to the genus Lasiopodomys (Abramson & Lissovsky, 2012). Morphological characteristics, such as karyotype (Gladkikh et al., 2016) and mating behavior (Zorenko & Atanasov, 2017), supported its current taxonomic status as L. gregalis. On the other hand, L. fuscus is nested in the genus Neodon Hodgson 1849 clade based on the longer length of ear and tail and greater number of inner angles in M1 and M3 compared with the genus Lasiopodomys (Liu et al., 2012a); moreover, CLOCK, BMA1, and Cytb gene sequences and their complete mitochondrial genomes supported this taxonomical status (Abramson et al., 2009a; Bannikova et al., 2010; Li, Lu & Wang, 2016a; Li et al., 2019; Liu et al., 2017). Recent studies have typically recognized Lasiopodomys as a separate genus that includes the species L. mandarinus and L. brandtii; L. gregalis was not widely accepted, whereas L. fuscus has been transferred to the genus Neodon and named Neodon fuscus.
According to fossils and molecular data, the genus Lasiopodomys originated and speciated during the Pleistocene epoch (∼2.58–0.012 million years ago (Ma)) when quaternary glaciations occurred in this period. Nuclear and mitochondrial phylogenetic estimates have shown that Lasiopodomys originated ∼2.4 Ma, whereas the division between L. gregalis and Lasiopodomys has been estimated to have occurred 1.8 Ma and that between L. mandarinus and L. brandtii was estimated at 0.5–0.95 Ma (Abramson et al., 2009b; Petrova et al., 2015; Li et al., 2017). However, chromosome analysis has shown that karyotype evolution has occurred between L. mandarinus and L. brandtii at ∼2.4 Ma, between Lasiopodomys and L. gregalis at 2.4 Ma, and between other Microtus species at 3 Ma (Gladkikh et al., 2016).
The species in the genus Lasiopodomys inhabit subterranean and aboveground environments and have recently become model species for comparative hypoxia adaptation (Dong et al., 2018; Sun et al., 2018). Species’ adaptation to low oxygen has been reported in numerous studies (Childress & Seibel, 1998; Dong et al., 2018; Nevo, 2013; Witt & Huerta-Sánchez, 2019), and most research has focused on animal models in an artificial environment or has compared them with subterranean rats to reveal the mechanisms of hypoxia (Ashur-Fabian et al., 2004; Malik et al., 2012; Malik et al., 2016). The differences in the environmental adaptability of proximal species are closely related to the historical events experienced during evolution, which play a key role in our understanding of the causes of current differences in life history among these species. However, the historical event that caused the Lasiopodomys species to adapt to a different environment has rarely been mentioned (Dong et al., 2018; Dong, Wang & Jiang, 2020).
Mitochondrial DNA are widely used to study the molecular ecology of animals because it is convenient and economical (Ballard & Rand, 2005; de Freitas et al., 2018; Kenechukwu, Li & An, 2018; Zhang et al., 2018). However, several studies have reported the limitations of mitochondrial DNA use (Galtier et al., 2009), such as recurrent horizontal transfer (Bergthorsson et al., 2003) and adaptive evolution (Bazin, Glémin & Galtier, 2006). The mitochondrial genome is involved in respiratory functions, which are closely associated with oxygen availability (Jain et al., 2016; Santore et al., 2002; Solaini et al., 2010).
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 high-throughput 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.
Materials and Methods
Material preparation and DNA sequencing
Total genomic DNA were extracted from the specimens of L. mandarinus (collected from 34°52′N, 113°85′E; Specimen No. LM023), L. brandtii (collected from 40°53′N, 116°38′E; Specimen No. LB003), and N. fuscus (collected from 34°9′N, 100°2′E; Specimen No. LF010) using the TIANamp Genomic DNA Extraction Kit (TIANGEN, DP304). All specimens were stored at the Animal Museum of Zhengzhou University. The Illumina NovaSeq 6000 (Illumina, San Diego, CA, USA) platform was used for sequencing the samples with a short-insert of 150 bp at ORI-GENE Company, Beijing (https://www.origene.com/).
Genome assembly and annotation
NOVOPlasty 3.6 was used for de novo assembly using the mitochondrial genome of L. mandarinus (GenBank no. JX014233) as a reference (Dierckxsens, Mardulyn & Smits, 2017). All mitochondrial genomes were annotated using GeSeq (Tillich et al., 2017), OGDRAW (Lohse et al., 2013), and GB2sequin (Lehwark & Greiner, 2019) in the MPI-MP CHLOROBOX integrated web tool (https://www.mpimp-golm.mpg.de/chlorobox), which contains the function of the HMMER package for protein-coding genes (PCGs) and ribosomal RNA (rRNA) (Finn, Clements & Edd, 2011), and tRNAscan-SE v2.0.3 for transfer RNAs (tRNAs) (Lowe & Eddy, 1997). Adenine–thymine (AT) skew was calculated as AT skew = (A − T) / (A + T), whereas guanine–cytosine (GC) skew was calculated as GC skew = (G − C)/(G + C). Circular maps were drawn using the CGView Server V 1.0 web tool (http://stothard.afns.ualberta.ca/cgview_server/) for L. mandarinus, L. brandtii, L. gregalis (GenBank no. MN199169), and N. fuscus (Grant & Stothard, 2008).
Molecular phylogenetic analysis and divergence time estimation
Phylogenetic analyses were performed on the whole mitochondrial genome sequences (Appendix S1). Besides the nine mitochondrial genomes that were acquired for the present study, five previously published mitochondrial genomes from L. mandarinus, L. gregalis, and N. fuscus were included; therefore, overall, 37 complete mitochondrial genome sequences from 23 species from the subfamily Arvivolinae were considered for phylogenetic analysis. Moreover, three species from Cricetulus Milne-Edwards 1867 were chosen as the outgroup. All these sequences were aligned using MAFFT v7.450 (Katoh & Standley, 2013). The nucleotide diversity of the PCGs of Lasiopodomys and Arvivolinae was determined using the DNASP v6.12.03 software (Rozas et al., 2017), and the best nucleotide substitution models were constructed using jMODELTEST 2.1.7 and selected using the Akaike information criterion (Darriba et al., 2012).
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).
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.
|Genes||Position (bp)||Strat/stop codon|
|L. brabdtii||L. mandarinus||L. gregalis||Neodon fuscus||L. brabdtii||L. mandarinus||L. gregalis||Neodon fuscus|
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).
|Species||contents||T||C||A||G||GC skew||AT skew|
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.781–2.58 Ma), with a possible most common ancestor of Lasiopodomys at ∼1.79 Ma (95% HPD values: ∼1.52–2.09 Ma). The split between L. brandtii and L. mandarinus was dated to the early Pleistocene period at ∼0.76 Ma (95% HPD values: ∼0.58–0.98 Ma), whereas the separation of both from L. gregalis was dated to the early Pleistocene epoch at 1.53 Ma (95% HPD values: ∼1.26–1.81 Ma). The estimated divergence event of N. fuscus and N. irene was found to be during the early Pleistocene epoch at 1.44 Ma (95% HPD Interval: ∼1.12–1.75 Ma).
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).
|Gene||Model||lnL||Models compared||Parameter Estimates||LRT ( P-value)|
|ATP6||Branch-model||A:One-ratio||−6321.793406||ω= 0.02625||p < 0.01|
|B:Omega = 1||−7802.578602||B vs A||ω=1|
|Branch-site model||Null||−6298.662308||null vs A||7 A 0.578||P<0.01|
|Model A||−6295.340396||60 I 0.987*|
|ATP8||Branch-model||A:One-ratio||−2156.664937||ω=0.16120||p < 0.01|
|B:Omega = 1||−2304.668614||B vs A||ω=1|
|Model A||−2092.811906||null vs A||NA|
|Cox1||Branch-model||A:One-ratio||−12806.04872||B vs A||ω=0.00534||p < 0.01|
|B:Omega = 1||−17316.74494||ω=1|
|Branch-site model||Null||−12712.24034||null vs A||57 I 0.779||0.077|
|Model A||−12710.67689||487 T 0.965*|
|Cox2||Branch-model||A:One-ratio||−5779.531238||ω=0.01386||p < 0.01|
|B:Omega = 1||−7512.338542||B vs A||ω=1|
|Model A||−5711.411586||null vs A||NA|
|Cox3||Branch-model||A:One-ratio||−6864.103188||ω=0.01989||p < 0.01|
|B:Omega = 1||−8741.926003||B vs A||ω=1|
|Model A||−6757.116417||null vs A||50 N 0.642||0.9065|
|62 V 0.517|
|203 F 0.593|
|B:Omega = 1||−12504.74705||B vs A||ω=1||p < 0.01|
|Model A||−10009.07827||null vs A||4 M 0.976*||0.1552|
|7 K 0.892|
|116 I 0.567|
|242 V 0.522|
|315 I 0.516|
|B:Omega = 1||−11391.13236||B vs A||ω=1||p < 0.01|
|Model A||−9015.745075||null vs A||NA||0.738|
|B:Omega = 1||−13190.51757||B vs A||ω=1||p < 0.01|
|Model A||−11268.21175||null vs A||11 F 0.747||1|
|14 F 0.816|
|31 I 0.845|
|95 T 0.837|
|122 I 0.856|
|207 I 0.845|
|220 H 0.867|
|228 K 0.847|
|235 N 0.860|
|241 L 0.858|
|B:Omega = 1||−4686.550566||B vs A||ω=1||p < 0.01|
|Model A||−3969.013478||null vs A||6 A 0.811||0.7962|
|14 S 0.790|
|20 V 0.861|
|108 Q 0.849|
|B:Omega = 1||−17886.34941||B vs A||ω=1||p < 0.01|
|Model A||−14856.89325||null vs A||NA||0.992|
|B:Omega = 1||−3775.223753||B vs A||ω=1||p < 0.01|
|Model A||−3151.8857||null vs A||NA||1|
|B:Omega = 1||−23436.47839||B vs A||ω=1||p < 0.01|
|Model A||−19737.31225||null vs A||194 E 0.512||0.9563|
|575 K 0.969*|
|B:Omega = 1||−5814.462234||B vs A||ω=1||p < 0.01|
|Model A||−4971.821283||null vs A||NA||1|
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) and Li 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, and Cytb (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 O2also 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.