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

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).

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 O₂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.