A new perspective on the taxonomy and systematics of Arvicolinae (Gray, 1821) and a new time-calibrated phylogeny for the clade

Background Arvicoline rodents are one of the most speciose and rapidly evolving mammalian lineages. Fossil arvicolines are also among the most common vertebrate fossils found in sites of Pliocene and Pleistocene age in Eurasia and North America. However, there is no taxonomically robust, well-supported, time-calibrated phylogeny for the group. Methods Here we present well-supported hypotheses of arvicoline rodent systematics using maximum likelihood and Bayesian inference of DNA sequences of two mitochondrial genes and three nuclear genes representing 146 (82% coverage) species and 100% of currently recognized arvicoline genera. We elucidate well-supported major clades, reviewed the relationships and taxonomy of many species and genera, and critically compared our resulting molecular phylogenetic hypotheses to previously published hypotheses. We also used five fossil calibrations to generate a time-calibrated phylogeny of Arvicolinae that permitted some reconciliation between paleontological and neontological data. Results Our results are largely congruent with previous molecular phylogenies, but we increased the support in many regions of the arvicoline tree that were previously poorly-sampled. Our sampling resulted in a better understanding of relationships within Clethrionomyini, the early-diverging position and close relationship of true lemmings (Lemmus and Myopus) and bog lemmings (Synaptomys), and provided support for recent taxonomic changes within Microtini. Our results indicate an origin of ∼6.4 Ma for crown arvicoline rodents. These results have major implications (e.g., diversification rates, paleobiogeography) for our confidence in the fossil record of arvicolines and their utility as biochronological tools in Eurasia and North America during the Quaternary.


INTRODUCTION
Arvicoline rodents (voles, lemmings, muskrats, and their extinct relatives) are the subject of a complex history of taxonomic and phylogenetic research (Conroy & Cook, 1999;Conroy & Cook, 2000;Galewski et al., 2006;Buzan et al., 2008;Robovský, ŘIčánková & Zrzavý, 2008;Fabre et al., 2012;Martínková & Moravec, 2012;Steppan & Schenk, 2017;Upham, Esselstyn & Jetz, 2019;Withnell, 2020;Abramson et al., 2021).Previous phylogenetic studies focused on subsets of Arvicolinae (i.e., North American Microtus; Conroy & Cook, 1999;Martínková & Moravec, 2012), used only one or two genes (Buzan et al., 2008), included arvicolines in larger studies investigating the phylogeny of all rodents (Fabre et al., 2012;Steppan & Schenk, 2017), or focused on mitochondrial genomes (Abramson et al., 2021).A comprehensive combined-evidence molecular and morphological dataset was generated by Robovský, ŘIčánková & Zrzavý (2008).Although a great deal of molecular data is available for arvicoline rodents, no study has synthesized all available molecular data to examine phylogenetic relationships across the group at both the generic and species levels-a new phylogeny has been warranted.Our objectives with this study are two-fold.(1) We wanted to consolidate a current and succinct understanding of the taxonomy and systematics of all the major clades of arvicolines through an analysis of the molecular data available on GenBank.(2) We wanted to use the fossil record to generate a time-calibrated phylogeny for the clade that can serve as the foundation that will allow for quantitative modeling of paleobiogeography and dispersal patterns of arvicolines across Beringia.To accomplish these objectives, we compiled the most taxonomically complete molecular dataset of global Arvicolinae to date (July, 2023), using published nuclear genes and mitochondrial markers, to provide a well-sampled molecular hypothesis of extant Arvicolinae.

Arvicoline taxonomy and systematics through time
For the sake of simplicity we have organized our summary of arvicoline taxonomy and systematics according to the Mammal Diversity Database (2023) of the American Society of Mammalogists.Different tribal affinities have been designated over the years and where relevant we have indicated this change.Tribal affinities in the headings reflect the current (August 2023) understanding recognized by the Mammal Diversity Database (2023).
Members of Lasiopodomys were considered by paleontologists to be the remnants of a group that was previously more speciose and widespread (Gromov & Polyakov, 1992;Repenning, 1992).We note that the fossil Lasiopodomys referred to by Repenning (1992) in North America is not the same as the extant Eurasian taxa, further adding to the taxonomic confusion of the genus (Repenning & Grady, 1988).Neontologists and paleontologists have recognized the morphological uniqueness of Lasiopodomys, but in one allozyme analysis, Lasiopodomys brandtii was grouped with Microtis fortis and Microtus ( =Stenocranius) gregalis, thus questioning the generic affinity of these species (Mezhzherin, Zykov & Morozov-Leonov, 1993).However, Musser & Carleton (2005) retained Lasiopodomys at the generic level, and Robovský, ŘIčánková & Zrzavý (2008) recognized Lasiopodomys as the sister to Stenocranius.Lemmiscus is a monotypic genus that was long considered a subgenus of Lagurus in order to segregate New World sagebrush voles from Old World steppe voles (Carroll & Genoways, 1980).Morphological and molecular data, however, indicate that Lemmiscus may be closely related to Microtus (Carleton, 1981;Modi, 1987;Abramson et al., 2021).
The taxonomy and systematics of Microtus are complicated and historically difficult to disentangle.Little consensus exists in the literature on how to treat generic-level identifications of fossil Microtus, partially because many hypotheses of Microtus relationships were based on tooth characters that have limited systematic potential and have undergone rapid evolutionary change (Guthrie & Matthews, 1971;von Koenigswald, 1980).Combined with the broad Holarctic distribution of the group, poor genetic sampling, and hypothesized recent origination and diversification, and the result has been taxonomic and systematic chaos.Currently, 60 species of Microtus are recognized within six subgenera; Blanfordimys, Euarvicola, Iberomys, Microtus, Pitymys, and Terricola (Kryštufek & Shenbrot, 2022;Mammal Diversity Database, 2023).Microtus (Blanfordimys) is a geographically isolated group of voles found in south-central Asia (e.g., Afghanistan) (Shenbrot & Krasnov, 2005).They have retained dental characters that have been interpreted as pleisiomorphic, but they have inflated auditory bullae and a mastoid region that is so enlarged that it almost projects beyond the occipital condyle, both of which have been interpreted as highly apomorphic (Gromov & Polyakov, 1992).This led some researchers to place them as a subgenus of Microtus (Gromov & Polyakov, 1992;Kryštufek & Shenbrot, 2022) while others gave them full generic distinction (Musser & Carleton, 1993).Bannikova, Lebedev & Golenishchev (2009) and Steppan & Schenk (2017) recovered Blanfordimys as sister to Euarvicola and this should be further examined.
Mictomicrotus is a recently named monotypic genus (Kryštufek & Shenbrot, 2022).Liu et al. (2007) originally placed M. liangshanensis in Proedromys (P.liangshanensis).However, morphology and phylogenetic studies demonstrated that liangshanensis is not the sister taxa to Proedromys bedfordi and therefore warrants a genus of its own (Chen et al., 2012;Steppan & Schenk, 2017).Voles of the genus Neodon are found throughout the mountainous regions of southern Asia.They have a long and complicated taxonomic history, but their close relationship to Microtus has been established, although systematic relationships of the genus relative to other arvicolines are still debated (Musser & Carleton, 2005;Pradhan et al., 2019).Recent taxonomic revision has seen the number of species belonging to Neodon grow (Liu et al., 2017;Pradhan et al., 2019).
Using morphological characters, Proedromys was hypothesized to be closely related to Microtus, but its diagnostic traits (massive cranium with wide, heavy, and grooved upper incisors and odd molars) were also used to support the hypothesis of a close relationship with extinct genera such as Allophaiomys (Gromov & Polyakov, 1992;Repenning, 1992).Molecular phylogenies also suggest that Proedromys is outside of and thus separate from Microtus, although with low support (Chen et al., 2012).Currently, Proedromys is considered a monotypic genus in the tribe Microtini and thought to be closely related to but separate from Microtus (Mammal Diversity Database, 2023).The narrow-headed voles, Stenocranius, have been considered as a subgenus of Microtus as well as Lasiopodomys (Kryštufek & Shenbrot, 2022).For many decades only one species (S. gregalis) was recognized, but recent molecular work has detected a cryptic species, S. raddei (Petrova et al., 2015;Petrova et al., 2016).Volemys consists of two high-altitude alpine species native to western Sichuan, China (Liu et al., 2017).Species of Volemys were previously placed in Microtus or were found to be closely related to Microtus, and previously published phylogenetic analyses of molecular data hinted that the distribution of Volemys may be relictual due to geographic (and correspondingly, genetic) isolation during the Late Pleistocene (Lawrence, 1982;Zagorodnyuk, 1990).

Tribe Pliophenacomyini (Arborimus and Phenacomys)
Similarities in dental morphology led some researchers to classify Arborimus as a subgenus of Phenacomys (Repenning & Grady, 1988), but others treated Arborimus as a separate genus (Musser & Carleton, 1993).Another study placed them together in the tribe Phenacomyini (Zagorodnyuk, 1990).Others placed Phenacomys with Phaiomys and other extinct genera (Repenning, Fejfar & Heinrich, 1990) or with the tribe Myodini (McKenna & Bell, 1997).Both Phenacomys and Arborimus have primitive molars that retain the plesiomorphic condition of retaining roots on molars, and they lack cementum in the reentrant angles on those molars; therefore some paleontologists argued that Phenacomys is an early relict lineage (Repenning, 1987).Currently there are two species of Phenacomys and three species of Arborimus recognized in two distinct but closely related genera united within the tribe Pliophenacomyini (Wilson, Lacher & Millermeier Jr, 2017).

Taxon sampling
Complete sampling of Arvicolinae has been historically challenging due to the high species diversity, global distribution of the clade, and the relative rarity of some species in museum collections.We attempted to sample all genera (n = 32) and species (n = 178) recognized by the Mammal Diversity Database ( 2023) in July 2023.We used Musser & Carleton (2005), Wilson, Lacher &Millermeier Jr (2017), andKryštufek &Shenbrot (2022) to inform our taxonomic coverage, however we used the Mammal Diversity Database ( 2023), the most current and widely accepted database of mammalian taxonomy, when calculating our taxonomic coverage.That resulted in three datasets; (1) a dataset of only taxa with mitochondrial data (n = 146), (2) a dataset of only taxa with nuclear data (n = 107), (3) a concatenated dataset that includes both mitochondrial and nuclear loci of n = 146 species of extant arvicolines, and is the most taxonomically complete dataset to date (August, 2023) for Arvicolinae (82% species and 100% generic coverage).Portions of this text were previously published as part of a thesis (https://repositories.lib.utexas.edu/bitstream/handle/2152/82563/WITHNELL-DISSERTATION-2020.pdf?isAllowed= y{&}sequence=1).
Sequences were aligned using the iterative refinement algorithm L-INS-I of MAFFT (Katoh & Standley, 2013).Aligned nexus files were imported into AliView (Larsson, 2014) and nuclear protein coding loci were checked for stop codons and trimmed where needed to ensure that they were in the proper reading frame for the first and third codon positions.PartitionFinder 2 (Lanfear et al., 2017) was used to partition the dataset (by codon position for the nuclear protein-coding genes) using the Akaike Information Criterion (AIC) (Burnham & Anderson, 2004).With all analyses we used GTR+ or GTR + + I molecular substitution models as suggested by PartitionFinder 2 (Lanfear et al., 2017) .For the specific model used with each partition see the results of the PartionFinder 2 analysis in Appendix S1B.

Phylogenetic analyses
We conducted Maximum Likelihood (ML) and Bayesian Inference (BI) analyses of the concatenated datasets, including mitochondrial only (Cytb and COI ), nuclear only (Ghr, IRBP/RBP3, BRCA1), and analyses of all five nuclear and mitochondrial markers, for a total of six phylogenetic analyses.We did not estimate gene trees for the nuclear exons or investigate a species tree approach because of uneven taxonomic sampling within each locus relative to the total taxonomic sample.With more comprehensive taxon sampling across loci, a species tree approach would be highly informative and could possibly help resolve some of the more challenging areas of the tree (e.g., within Microtus).

Time-calibrated analyses
We conducted a time-calibrated BI analysis in MrBayes 3.2.7 using our concatenated mitochondrial and nuclear datasets.We used a birth-death model and an independent gamma rate relaxed-clock (igr), where each branch has an independent rate drawn from a gamma distribution that was empirically derived in MrBayes.The MCMC chain was run for 3.0 × 10 7 generations (sampled every 1,000 generations) for two runs each with four chains.A temperature of 0.1 was implemented and the first 30% of the data were discarded as burn-in.Results of the analyses were visualized in Tracer v1.7 (Rambaut et al., 2018) to ensure runs had reached stationarity and that the effective sample size (ESS) was >200 for all model parameters.
Microtus is one of the most diverse and rapidly evolving mammalian genera (Triant & De Woody, 2006).Many phenotypic characters are convergent among distantly related species, and high genetic variation has been attributed to karyotypic differentiation, with diploid chromosomal numbers ranging from 17 to 64 (Triant & De Woody, 2006).Triant & De Woody (2006) documented that Microtus sensu stricto has a time-corrected mitochondrial rate of nucleotide substitution of 0.08 substitutions per site per million years.This is higher than most other mammals (e.g., Pan, Bos, Ursus) and obviously would affect divergence time analyses (Triant & De Woody, 2006).We therefore calculated the rate of evolution for the mitochondrial only, nuclear only, and concatenated mitochondrial and nuclear datasets using R code provided in Gunnell et al. (2018).This code uses a user inputed tree-age along with branch lengths to calculate the rate of evolution within the tree.The mitochondrial rate of evolution was 0.09 substitutions per site per million years, the nuclear rate of evolution was 0.01 substitutions per site per million years, and the concatenated dataset was 0.065 subtitutions per site per million years.This code is included in our Dryad submission (Appendix S1C).Since we used the concatenated mitochondrial and nuclear dataset for time-calibration we chose to use the substitution rate of 0.065 substitutions per site per million years as the mean clock rate prior.We used the MrBayes command 'prset clockratepr' with a mean of −2.72 (natural log of 0.065) and a standard deviation of 0.12 as calculated in the code from Gunnell et al. (2018).

Node calibration selection
We used four internal node calibrations and a root calibration in our divergence time analyses.Calibrated nodes were constrained as monophyletic and these nodes were selected after non-calibrated phylogenies were produced.For all nodes, there were no suitable fossils available to help establish calibration maxima, so we used exponential calibration priors for each node.For each calibration, the fossil age was used as the offset.R scripts for calculating a suitable mean are in Appendix S1C.

Calibration 1: Cricetidae (outgroup) (tree root)
We chose as outgroups three muroid rodents previously found to be closely related to Arvicolinae (Fabre et al., 2012;D'Elía, Fabre & Lessa, 2019).These three species belong to the subfamilies Cricetinae (Cricetus cricetus, Mesocricetus auratus) and Neotominae (Neotoma fuscipes).The split between Arvicolinae and Cricetinae is reported to have occurred during the middle Miocene (Fabre et al., 2012).The split between (Arvicolinae, Cricetinae) and Neotominae was hypothesized to be during the early-Miocene (Fabre et al., 2012;Steppan & Schenk, 2017).We chose to calibrate the root of our tree using the oldest purported fossil neotomid, Lindsaymys sp.cf.L. takeuchii (Kelly & Whistler, 2014;Martin & Zakrzewski, 2019).This calibration is anchored by a lower first molar (m1) housed at the Los Angeles County Museum (LACM 157168).LACM 157168 is diagnosed as belonging to Lindsaymys based on.(1) a moderately hyposodont molar that is smaller than the more temporally recent specimens (2) an m1 with an anteroconid that is not bifurcated and positioned close to the metaconid (3) a metalophulid that connects to the protolophulid I at the junction with the anterolophid (4) presence of an entoconid spur and (5) a moderately deep valley between the metaconid and the lingual edge of the anteroconid and anterlophid (Kelly & Whistler, 2014).This last feature is key to distinguishing it from other contemporary taxa such as Abelmoschomys, Antecalomys, Prosigmodon, Bensonomys, Baiomys, Symmetrodontomys, and Jacobsomys (Kelly & Whistler, 2014).LACM 157168 was found in locality LACM 5720 which is thought to be Latest Clarendonian to early Hemphillian (C13-Hh1, ∼9.2−8.7 Ma) (Kelly & Whistler, 2014).LACM Locality 5720 is found in the El Paso Mountains within the Dove Spring Formation of the western Mojave Desert, California (Kelly & Whistler, 2014).The site lies underneath Dove Spring Ash number 16 dated via Ar/Ar at 8.5 ± 0.13 Ma and above Dove Spring Ash number 15 dated via fission tract at 8.4 ± 1.8 Ma (Whistler et al., 2009).Biochronology and paleomagnetics place this site in Chron C4A, with a maximum age of ∼9.2 Ma (Kelly & Whistler, 2014).We chose to use 8.8 Ma as a conservative estimate of the minimum age of the site, since it is a middle point between 8.5 Ma and 9.2 Ma (Whistler et al., 2009;Kelly & Whistler, 2014).The soft maximum for the node was based on the divergence time analysis of Steppan & Schenk (2017).We used an offset exponential distribution with a minimum age of 8.8 Ma and a mean of 10.8 Ma. 10.8 Ma was chosen as the mean because it produced an 95% upper bound of the distribution at 15 Ma.

Calibration 2: Lemmini node
The earliest North American bog lemmings (Synaptomys) are from the Hagerman Fossil Beds National Monument, Idaho (Mictomys = Synaptomys vetus; (Ruez Jr & Gensler, 2008).The offset for this node is anchored by a right m1 (lower first molar) housed at the Idaho Museum of Natural History (IMNH 67002/39517) that has radiometric age control (Ar-Ar) of a basaltic tephra located 30 m above the site and dated at 3.79 ± 0.03 Ma (Hart & Brueseke, 1999).Interpolation of depositional rates indicates that the age of the fossil from IMNH locality 67002 is ∼3.95 Ma (Hart & Brueseke, 1999).IMNH 67002/39517 was identified as Mictomys vetus by having evergrowing molars with cementum in the reentrant angles.The m1 also has a posterior loop with three triangles, and an anterior loop (Ruez Jr & Gensler, 2008).Triangles 1 and 2 are broadly confluent with the anterior loop and triangle three is joined by the anterior loop near the midline.Triangles 1 and 3 are nearly twice the width of triangle 2. Because Synaptomys was paraphyletic in some of our uncalibrated analyses, we used this fossil to calibrate the crown lemming node instead.We used an offset exponential distribution with a minimum age of 3.95 Ma and a mean of 4.74 Ma (see Appendix S1C for this calculation).

Calibration 3: Ondatrini node
The oldest known species of this clade, Ondatra minor, is found in the Hagerman Formation in Hagerman, Idaho at ∼3.6 Ma (Hibbard, 1959).All of the fossils at Hagerman are constrained between two lava flows and ash units that have yielded ages of 4.0 Ma to 3.2 Ma using Ar-Ar dating methods (McDonald, Link & Lee, 1996).We anchored the Ondatra + Neofiber node using a left m1 tooth of Ondatra minor (USNM 21830) from Hagerman.This m1 was identified as Ondatra minor by its relatively large size as well as being rooted and having a posterior loop, five alternating triangles, with a fifth triangle opening broadly into the anterior loop (Hibbard, 1959).We used an offset exponential distribution with a minimum age of 3.2 Ma and a mean of 4.9 Ma.The age of 3.2 Ma was chosen because it is the most conservative estimate of the age of the two ash layers described from Hagerman Idaho and deposition interpolation information were not available for the locality.

Calibration 4: Pliophenacomyini node
Extant voles of the genera Phenacomys and Arborimus are today restricted to North America.Eurasian specimens of Phenacomys were identified from Krestovka, Kolyma Lowland Russia (Sher et al., 1979;Zazhigin, 1997) and Romanovo 1c, Western Siberia, Russia (Smirnov, Bolshakov & Borodin, 1986;Borodin, 2012).Recently, a new species (Phenacomys europaeus) was described from Europe in Zuurland, the Netherlands, and dated at ∼2.1 Ma based on biochronology (van Kolfschoten, Tesakov & Bell, 2018).The oldest known record of Phenacomys, P. gryci, (type locality in the Gubik Formation) is from the Fish Creek fauna of Alaska.The Fish Creek Fauna is in the Gubik Formation, which is an alternating marine and coastal plain sedimentary unit.The Fish creek Fauna is dated at ∼2.4 Ma using amino acid racemization ratios, a reversed polarity zone, and the presence of the ancestral sea otter Enhydrion and the arvicoline rodent Plioctomys mimomiformis (Carter et al., 1986;Repenning et al., 1987;Repenning & Brouwers, 1992).We calibrated the (Phenacomys + Arboriumus) node based on the type specimen of Phenacomys gryci (a left m1 housed at the United States National Museum USNM 26495).This fossil was assigned to Phenacomys gryci by having a rooted m1 that lacked cementum in the reentrant angles.It also possesses a posterior loop, five asymmetrical alternating triangles with a ''Mimomys Kante'' on triangle four, and a complex anterior loop (Repenning et al., 1987).This node was calibrated using an offset exponential distribution with a minimum age of 2.4 Ma and a mean of 3.27 Ma (see Appendix S1C for this calculation).

Calibration 5: Ellobiusini node
The timing of the origination and diversification of Microtus and its close relatives has been repeatedly contested among paleontologists (e.g., Repenning, 1992;Martin & Tesakov, 1998).It was argued that the genus Allophaiomys gave rise via anagenetic evolution to what is recognized today as Microtus (Martin & Tesakov, 1998), but that hypothesis is controversial (e.g., Bell et al., 2004;Bell & Bever, 2006).The oldest Allophaiomys with external age control is from Hansen Bluff (Colorado) and dated at 1.9 Ma (Rogers et al., 1992).The earliest occurrence of Microtus, as defined by Repenning (1992), was long thought to be from the Anza-Borrego Desert of California (Zakrzewski, 1972), possibly from 1.4 to 1.6 Ma (lacking firm age control).Unfortunately, the specimens from Anza-Borrego had questionable field data; one specimen was found in a fault block and from a different area in the park than originally reported, and a second specimen could not definitively be assigned to Microtus (Bell & Bever, 2006;Murray, Ruez Jr & Bell, 2011).The oldest known Microtus is, therefore, found in the type Irvington Fauna from California dated to 1.21 Ma based on paleomagnetic data (Bell & Bever, 2006).
Fossil evidence from Ellobius, an early diverging member of the clade that includes Microtus, was used here to calibrate the node.The oldest fossils of Ellobius are from the Late Pliocene of Kazakhstan and Tajikistan (Lytchev & Savinov, 1974;Zazhigin, 1988) and the Northern Caucasus (Tesakov, 2004).We chose to use a fossil mandible with m1-m3 (Paleontological Institute, Russian Academy of Sciences M-2049/58-KB) of Ellobius primigenius from Central Asia (Lytchev & Savinov, 1974).This fossil possesses rooted teeth with relatively high crowns, a posterior loop, five alternating triangles, and an anterior loop consistent with Ellobius (Lytchev & Savinov, 1974).This mandible is part of the Kiikbai Fauna of Kazakhstan dated using biochronology (the occurrence of Hypolagus brachygnathus, Ochotonoides complicidens, and Mimomys pliocaenicus) to the Pliocene at ∼2.4 Ma in the Matuyama Chron (Sotnikova, Dodonov & Pen'Kov, 1997).The Kiikbai Fauna is described from the southern flank of the Ilian depression in the Alatau mountains and placed in the European middle Villafranchian land mammal age (Sotnikova, Dodonov & Pen'Kov, 1997).We used an offset exponential distribution with a minimum age of 2.4 Ma and a mean of 3.27 Ma (see Appendix S1C for this calculation).

Non-clock analyses
Six phylogenetic analyses were conducted on the concatenated dataset that included either all or a subset of the 147 sequences from Cytb; 64 from COI ; 90 from Ghr; 105 from IRBP/RBP3; and 69 from BRCA1.110 species (74%) included both mitochondrial and nuclear data.39 species (26%) had only mitochondrial data.GenBank accession numbers are in Appendix S1A.Results from Maximum Likelihood and Bayesian analyses were similar or identical except where discussed below.

Bayesian inference (BI) results
A tribe-level summary of the BI tree including only the mitochondrial loci (Cytb and COI) with posterior probability values is presented in Fig. 1B.High support (>95 PP) was found for 11 (32%) of the nodes (Fig. 1B).A species-level tree is presented in Appendix S1G.High support was found for 80 (54%) of the nodes (Appendix S1G).A tribe-level summary of the BI tree including only the nuclear loci (Ghr, IRBP/RBP3, and BRCA1) with posterior probability values is presented in Fig. 2B.High support was found for 20 (74%) of the nodes (Fig. 2B).A species-level tree is presented in Appendix S1H.High support was found for 68 (63%) of the nodes (Appendix S1H).Finally a tribe-level summary of the BI tree including all five loci (Cytb, COI, Ghr, IRBP/RBP3, and BRCA1) with posterior probability values is presented in Fig. 3B.High support was found for 17 (57%) of the nodes (Fig. 3B).A species-level tree is presented in Appendix S1I.High support was found for 88 (60%) of the nodes (Appendix S1I).

Major clade systematic results
The tribe Arvicolini was variably placed in the different analyses (Figs.1-3).In the mitochondrial analyses (Fig. 1) Arvicolini is weakly placed near the root of the arvicoline tree.In the nuclear and combined analyses (Figs.2-3) it is placed with stronger support in the derived clade that includes Ellobiusini, Pliomyini, Lagurini, and Microtini.Clethrionomyini is consistently found in all analyses to be a large clade in the middle of the arvicoline tree (Figs.1-3).Within this tribe there are two clades that are found in all analyses.(1) A clade that consists of Caryomys, Eothenomys, and Anteliomys (Figs. 1-3), and (2) A clade that consists of Craseomys, Clethrionomys, and Alticola (Figs. 1-3).Anteliomys is paraphyletic with strong support with respect to Eothenomys in the nuclear analysis (Fig. 2) but monophyletic and the sister to Eothenomys in the other analyses (Figs. 1 and 3).The paraphyly/monophyly of Clethrionomys and Alticola vary among analyses (Figs.1-3).Clethrionomys is paraphyletic in the ML and BI mitochondrial and nuclear analyses as well as the combined BI analysis (Figs. 1-3).In the combined ML analysis, Clethrionomys was monophyletic (Fig. 3A).Alticola was found to be paraphyletic in the mitochondrial analyses (Fig. 1) and the combined ML analysis (Fig. 3A) and monophyletic in the nuclear analysis (Fig. 2) and combined BI (Fig. 3B).
Dicrostonychini is weakly inferred as the sister to Ondatrini in the mitochondrial analyses (Fig. 1), and strongly supported as the sister to Pliophenacomyini in the nuclear and combined analyses (Figs.2-3).Ellobiusini is consistently found at the base of the large clade that includes Lagurini, Arvicolini (in the nuclear and combined analyses), Pliomyini (in the nuclear and combined analyses), and Microtini (Figs. 1-3).Lagurini is found to be the sister to Microtini in the mitochondrial analyses (Fig. 1) but sister to the  -3), but relationships within the tribe vary across analyses.Myopus is the first diverging member of Lemmini in the mitochondrial analysis (Fig. 1 and Appendices S1D and S1G), but earliest diverging 'true-lemming' in the nuclear and combined analyses (Figs.2-3 and Appendices S1E, S1H and S1F, S1I).Synaptomys is paraphyletic or polyphyletic in the mitochondrial analyses (Fig. 1) and combined ML analysis (Fig. 3A) but monophyletic in the nuclear analyses (Fig. 2) and combined BI analysis (Fig. 3B).
Microtini was consistently the most crownward tribe of arvicolines (Figs.1-3).Chionomys is sister to the rest of Microtini in the nuclear and combined analyses (Figs.2-3), but in the mitochondrial analyses (Fig. 1), a split between Microtus ochrogaster and the rest of Microtini was the first divergence.Volemys are early diverging within Microtini, and are polyphyletic in the ML mitochondrial analysis (Fig. 1A), nuclear analyses (Fig. 2), and combined analyses (Fig. 3).However, Volemys was found to be monophyletic in BI mitochondrial analysis (Fig. 1B).In all analyses, Proedromys was placed as sister to Volemys musseri (Figs.1-3).Lasiopodomys is the sister to Stenocranius in the BI mitochondrial only (Fig. 1B) and the nuclear and combined analyses (Figs.2-3).In the ML mitochondrial only analysis, Lasiopodomys is sister to Alexandromys, albeit weakly (Fig. 1A).In the ML mitochondrial only analysis, Stenocranius was sister to a paraphyletic Neodon (containg a nested Volemys millicens) (Fig. 1A), but in the others it is sister to Lasiopodomys (Fig. 1B-3).Neodon is found to be monophyletic in the mitochondrial BI analysis (Fig. 1B), nuclear ML (Fig. 2A) and paraphyletic in the other analyses (Figs.1A,  2B, and 3).Mictomicrotus is deeply nested within Microtini (Figs. 1A, 2 and 3) or is nested within Microtus (Fig. 1B).Alexandromys is one of the most nested members of Microtini (Figs. 1-3).Microtus is the most nested arvicoline and was monophyletic in the nuclear and combined analyses (Figs.2-3) but paraphyletic or polyphyletic in the mitochondrial analyses (Fig. 1).Endemic species of North American Microtus are monophyletic in the nuclear and combined analyses (Figs.2-3).
Ondatrini is a rootward tribe of arvicolines in all analyses, but its sister relationships vary (Figs.1-3).In the mitochondrial analyses Ondatrini is sister to Dicrostonychini (Fig. 1), and in the nuclear and combined analyses it is sister to Pliophenacomyini + Dicrostonychini + Lemmini (Figs. 2-3).Pliophenacomyini is recovered as an early diverging arvicoline sister to Dicrostonychini in the nuclear and combined analyses (Figs.2-3).Pliomyini is found towards the base of the arvicoline tree in the mitochondrial analyses (Fig. 1), but within the clade that includes Ellobiusini, Lagurini, Arvicolini, and Microtini in the nuclear and combined analyses (Figs.2-3).Prometheomyini is found at the base of the arvicoline tree in the nuclear and combined analysis (Figs. 2-3), but at the base of the clade that includes Ellobiusini, Lagurini, and Microtini in the mitochondrial analyses (Fig. 1).

Time-calibrated analysis
A tribal-level time-calibrated majority rule consensus tree of the combined dataset is presented in Fig. 4. A species-level tree is presented in Appendix S1J.105 nodes (72%) had posterior probability values >0.95.

Divergence-time results
Crown arvicoline rodents were inferred to have diverged ∼6.4 Ma.For a list of all the major clades and their divergence estimates see Table 1.
There are some topological differences based on marker set inclusion (e.g., nuclear, mitochondrial, or combined).Below we outline the implications of our results for the taxonomy and evolutionary understanding of Arvicolinae.
Prometheomyini was deeply nested (though with weak support) in the mitochondrialonly analyses (Fig. 1), and sister to arvicolines besides Lemmiscus in the combined ML analysis.Therefore, we stress caution in considering Prometheomys as the most basal arvicoline.Discrepancy between the mitochondrial and nuclear loci warrants further investigation (Pardiñas et al., 2017;Kryštufek & Shenbrot, 2022).

Systematic relationships within Lemmini
Bog lemmings (Synaptomys) and 'true lemmings' (Myopus, Lemmus) are consistently placed in a clade at or near the base of the arvicoline tree.One previous study placed lemmings as sister to all other arvicolines (Abramson et al., 2021), whereas most have found lemmings to be near but not at the base of the tree (Galewski et al., 2006;Buzan et al., 2008;Robovský, ŘIčánková & Zrzavý, 2008;Fabre et al., 2012;Steppan & Schenk, 2017;Upham, Esselstyn & Jetz, 2019).Studies that sampled both species of extant bog lemming found that Synaptomys (as defined by Musser & Carleton (2005)) is paraphyletic with respect to 'true lemmings' (Buzan et al., 2008;Fabre et al., 2012;Steppan & Schenk, 2017;Upham, Esselstyn & Jetz, 2019).Our results unfortunately do not add any clarity on whether bog lemmings should be one or two genera, because 50% of our trees recover a weakly supported clade.However, it should be noted that analyses based either entirely or predominantly on nuclear markers recovered a weakly supported clade, while those that included mitochondrial-only data tended not to.Based on dental morphology, many paleontologists have considered the northern bog lemming, Synaptomys borealis (following Musser & Carleton, 2005), to be a member of a distinct genus known as Mictomys (Fejfar & Repenning, 1998;Repenning & Grady, 1988).Musser & Carleton (2005) and Pardiñas et al. (2017) argued that there may be enough evidence to place Mictomys and Synaptomys as distinct genera, but they tentatively kept them in the same genus.
Our nuclear and combined mitochondrial and nuclear analyses are in agreement with Abramson et al. (2021) that Dicrostonyx, Phenacomys, and Arborimus, are united in a clade that Abramson et al. (2021) called Dicrostonychini.However, we have followed the tribal terminology of Kryštufek & Shenbrot (2022) in which Dicrostonychini (Dicrostonyx only) and Pliophenacomyini (Phenacomys + Arborimus) are closely related but separate tribes.
Lemmings and ((Phenacomys, Arborimus), Dicrostonyx) were consistently placed in a clade at or near the base of Arvicolinae in our analyses.Although it is true that Dicrostonyx are not sister to Lemmini, they are part of a more inclusive clade that is often plaed as sister to Lemmini.

Systematic status of Hyperacrius
The systematic position of the Subalpine Kashmir Vole (Hyperacrius fertilis) has been relatively understudied.Using morphology alone, Hyperacrius was hypothesized to be closely related to Alticola (Hinton, 1926), or as a member of the tribe Clethrionomyini (Gromov & Polyakov, 1977).Kohli et al. (2014) included Hyperacrius for the first time in a molecular analysis, and its relationship to Clethrionomyini was doubted.We found weakly-supported close relationships between Hyperacrius and Clethrionomyini (mitochondrial only and combined ML) and near Pliomyini (combined BI).Recently, Hyperacrius was hypothesized to be the earliest diverging member of what was previously considered Arvicolini (i.e., Microtini minus Arvicola) (Abramson et al., 2020;Abramson et al., 2021).Given that only cytb data were available for Hyperacrius our results should be considered equivocal and more sampling is needed.

Systematic status of Lemmiscus curtatus
Until recently, the sagebrush vole, Lemmiscus curtatus, had not been included in molecular phylogenies (Steppan & Schenk, 2017;Abramson et al., 2021).Those two analyses produced conflicting results for the systematic position of Lemmiscus.Both studies found that Lemmiscus and Microtus are not sister taxa.Steppan & Schenk (2017) placed Lemmiscus as sister to Arvicola, whereas Abramson et al. (2021) placed it as sister to Chionomys.Some of our results, though weakly supported, are similar to Steppan & Schenk (2017), probably because we used the same genetic data for Lemmiscus for our phylogeny, whereas Abramson et al. (2021) used an entire mitochondrial genome of Lemmiscus.Our combined ML analysis placed Lemmiscus as sister to all other arvicolines whereas our other analyses place it close to Arvicolini.More data from transcriptomes or nuclear genes will help to refine the systematic position of this species, so we do not make any taxonomic recommendations based on the systematic position of Lemmiscus in our analyses.
Arvicolini has been used to describe a large number of genera including Arvicola, Chionomys, Neodon, Lasiopodomys, and Microtus.We recognize the tribe Arvicolini as having a single genus, Arvicola, based on the Mammal Diversity Database (2023).We inferred a monophyletic Arvicolini with a weak sister taxon relationship to Lemmiscus curtatus (in mitochondrial only BI), similar to the results of Steppan & Schenk (2017).That is an interesting biogeographic result given the large distance between the extant members of these genera.We also recovered Arvicolini as the sister clade to Microtini (in both nuclear and combined analyses) and as in Abramson et al. (2021).
Lagurini was placed as the sister to Pliomyini (Steppan & Schenk, 2017), as the sister to Arvicola (Abramson et al., 2021), and as sister to Microtini (Fabre et al., 2012).Our nuclear only and combined results are most similar to Fabre et al. (2012) and Abramson et al. (2021) in that we found Lagurini to be closely related to Arvicolini and Microtini.Pliomyini was found to be sister to Lagurini (Steppan & Schenk, 2017), and as the sister to Arvicola (Fabre et al., 2012;Abramson et al., 2021).Regardless of its sister taxon relationship all three of these studies recovered Pliomyini close to Ellobiusini, Arvicolini, and Lagurini.Our analyses placed Pliomyini as (1) sister to Arvicolini and Clethrionomyini wth mitochondrial data only, or (2) as closely related to Arvicolini, Lagurini, Ellobiusini, and Microtini in the nuclear only and combined analyses.Based on those results, Pliomyini is part of the major radiation that includes Ellobiusini, Lagurini, Arvicolini, and Microtini.

Systematic status of Microtini
Microtini is the largest tribe of arvicoline rodents.Recently, the Mammal Diversity Database ( 2023) recognized Microtini as distinct from Arvicolini based on Abramson et al. (2021), who found that Arvicola was not in a clade with the other members of Microtini.We found the same result (see Figs. 1-3).Therfore we retain the genera Chionomys, Proedromys, Volemys, Mictomicrotus, Neodon, Lasiopodomys, Stenocranius, Alexandromys, and Microtus within a monophyletic Microtini.
Historically, Chionomys was placed in Arvicola, in Microtus, its own genus, or as a subgenus of Microtus (Yannic et al., 2012).This complicated history can be attributed in part to the fragmented geographic distribution and isolation in high alpine environments of Chionomys.Jaarola et al. (2004) used analyses of cytb to solidify Chionomys as a valid genus separate from Microtus, and several other studies placed Chionomys as a nested member of what is now considered Microtini, but outside of Microtus (Galewski et al., 2006;Robovský, ŘIčánková & Zrzavý, 2008;Fabre et al., 2012;Abramson et al., 2021).Our results further support that hypothesis and provide evidence to suggest that Chionomys is a basal member of Microtini.Proedromys has historically been thought to be closely related, or even included in, Microtus (Ellerman & Morrison-Scott, 1951;Gromov & Polyakov, 1977;Musser & Carleton, 2005).We concur with other phylogenetic studies (Fabre et al., 2012;Steppan & Schenk, 2017;Abramson et al., 2021) that place Proedromys in Microtini and close to but outside of Microtus.We consistently found Volemys musseri as sister to Proedromys with strong support (see Figs. 1-3).Given the geographical overlap of the two genera that result is not suprising and recapitulates previous studies (Steppan & Schenk, 2017;Upham, Esselstyn & Jetz, 2019).Mictomicrotus is a recently named genus that is monotypic (M.liangshanensis) and was previously included in Proedromys (Liu et al., 2007;Kryštufek & Shenbrot, 2022;Steppan & Schenk, 2017).Our results further support the hypothesis that Mictomicrotus is a distinct genus.
Recent systematic and taxonomic work altered our understanding of Microtus (Bannikova et al., 2010;Liu et al., 2012;Liu et al., 2017;Pradhan et al., 2019).For example, Neodon was previously recognized as its own genus or as a subgenus of Microtus (see Musser & Carleton (2005); Pardiñas et al. (2017)), or was placed in Pitymys (Ellerman & Morrison-Scott, 1951).It has become increasingly apparent that Microtus was historically used as a taxonomic garbage bin.Several species previously allocated to Microtus are now placed in Neodon (N.leucurus, N. clarkei, and N. fuscus) (Pradhan et al., 2019;Abramson et al., 2021).We found a paraphyletic Neodon with respect to Volemys millicens in most analyses (mitochondrial only ML, nuclear only BI, and combined ML and BI), and a monophyletic Neodon that is sister to Alexandromys in the others (mitochondrial only BI and nuclear only ML).The relationship between Volemys millicens and Neodon warrants further exploration.
Volemys was considered a distinct genus or a subgenus of Microtus, with Musser & Carleton (2005) recognizing two species, V. musseri and V. millicens.That classification is based on morphology alone, and the monophyly of Volemys has not been supported by molecular datasets (Jaarola et al., 2004;Steppan & Schenk, 2017;Upham, Esselstyn & Jetz, 2019).We inferred a weakly monophyletic Volemys in one analysis (mitochondrial only BI).We suggest further systematic study of Volemys and further examination to understand potential morphological homoplasy between V. millicens and V. musseri.
Voles placed within Lasiopodomys and Stenocranius have a long history of taxonomic change, with some studies placing them within Microtus, while others have placed them all within Lasiopodomys (see Musser & Carleton, 2005;Pardiñas et al., 2017;Kryštufek & Shenbrot, 2022).We follow the taxonomy recognized by Kryštufek & Shenbrot (2022) and the Mammal Diversity Database (2023) in which there are two species of Lasiopodomys and two species of Stenocranius.All but one of our analyses (mitochondrial only ML) found Lasiopodomys and Stenocranius to be sister taxa with medium to strong support, in agreement with past studies (Fabre et al., 2012;Steppan & Schenk, 2017;Abramson et al., 2021).
Alexandromys was used for holartic grass voles for more than a century, but its usage has waxed and waned (Kryštufek & Shenbrot, 2022).Voles that the Mammal Diversity Database (2023) assign to Alexandromys have historically been placed in Microtus, Iberomys (Gromov & Polyakov, 1977), Neodon (Musser & Carleton, 2005), and in various subgeneric groups within Microtus (Kryštufek & Shenbrot, 2022).Phylogenetic analyses, however indicate the validity of Alexandromys as a genus separate from Microtus.Steppan & Schenk (2017) and Abramson et al. (2021) found a monophyletic Alexandromys but Fabre et al. (2012) found the genus to be paraphyletic with respect to Volemys musseri.All of our analyses used a robust sampling within Alexandromys and all found a monophyletic genus that is closely related to but outside of crown Microtus (see Figs. 1-3).Given the aforementioned issues with Volemys it is likely that Alexandromys is a monophyletic genus whose recognition helps to make Microtus monophyletic.
Voles of the genus Microtus are frequently studied but have presented a long-term systematic enigma (e.g., Conroy & Cook, 2000;Fabre et al., 2012;Abramson et al., 2021).Microtus is one of the most rapidly evolving lineages of rodents and contains over sixty extant species (Mammal Diversity Database, 2023).To attempt to clarify the taxonomy of Microtus, researchers have used subgenera such as Pedomys, Alexandromys, Terricola, Iberomys, Agricola, and Neodon, but these subgenera are variably considered genera by different authors.Thus, the genus Microtus was in need of being redefined (Barbosa et al., 2018;Abramson et al., 2021).Recent taxonomic work concerning Alexandromys, Neodon, Lasiopodomys, Stenocranius, Chionomys, Proedromys, and Mictomicrotus has helped to address this issue and our results (a monophyletic in our nuclear only and combined analyses; see Figs. 2-3) further help to support those changes (Kryštufek & Shenbrot, 2022).This is a key step to clarifying this complicated genus and getting closer to being able to have an answer for 'What is Microtus?'

Systematic status of North American Microtus
The species of Microtus that are endemic to North America have in the past been inferred as a clade (Conroy & Cook, 1999;Conroy & Cook, 2000;Upham, Esselstyn & Jetz, 2019;Abramson et al., 2021).Our mitochondrial only results conflict with that hypothesis (see species level trees in Appendices S1D and S1G) in that North American Microtus is paraphyletic with respect to M. cabrerae (a species endemic to the Iberian Penninsula), while our nuclear only and combined analyses (see species level trees in Appendices S1E, S1H and S1F, S1I) found a weakly-supported clade.The relationship between M. cabrerae and the North American Microtus should be further examined using transcriptome and genomic data.Whether or not North American Microtus is monophyletic has important implications for paleobiogeography and our understanding of the fossil record in relation to North American Land Mammal Ages (NALMAs), which is largely based on the immigration of voles and other mammals into North America at various times during the Pliocene and Pleistocene (Bell et al., 2004).

Diversification of arvicolines
Time-calibrated phylogenies rely on several important factors to assure their accuracy and reproducibility.Node calibrations must leverage appropriate fossils to avoid erroneous divergence times.As recommended by Parham et al. (2012), care must be taken to insure that a fossil that is used for a calibration has (1) a museum number, (2) a systematic position established with an apomorphy based diagnosis or phylogenetic analysis, (3) reconciliation of morphological and molecular data with respect to the fossil's position, (4) detailed locality and geological data, and (5) reference to a published radiographic age.Following this protocol insures that a calibration is as accurate as possible, given the understanding of the fossil record at the time of the analysis.Studies have also demonstrated that secondary calibrations can dramatically affect the age of a node (Schenk, 2016).Secondary calibrations can certainly be used to help inform node calibrations, but we advise the use of primary fossil calibrations, when possible.In our divergence-time analysis we chose to implement a node-dating approach rather than tip-dating or the fossilized birth-death model (FBD) (Heath, Huelsenbeck & Stadler, 2014), because there are currently so few phylogenetically constrained arvicoline fossils.We acknowledge that the use of single fossils for node age minima and maxima (node-dating) is potentially problematic (Heath, Huelsenbeck & Stadler, 2014).In the future when there are more phylogenetically constrained fossils, we recommend a reexamination of this issue using FBD methods.
We argue that one of the strengths of our paper is our justification of our fossil calibrations following Parham et al. (2012).Our results are largely congruent with the published fossil record and our divergence time estimates are younger (closer to the known fossil record) than those produced by Steppan & Schenk (2017), Upham, &Jetz (2019), andAbramson et al. (2021).This is not surprising given the density of fossil calibrations that we used across the arvicoline tree.We estimated a median age of crown Arvicolinae at ∼6.4 Ma, which is slightly younger than the age (7.4 Ma) inferred by Abramson et al. (2021).These ages are very close to the age of the earliest known possible arvicoline rodent fossils (∼9-8 Ma, which are probably Pliocene in age (e.g., Pannonicola, Microtoscoptes, or Goniodontomys; Maul et al., 2017).Abramson et al. (2021) included a calibration at the crown Arvicolinae node (7 Ma FAD, 10 Ma max age) whereas we did not, potentially accounting for some discrepancy in our results.We did not use a calibration at the root of Arvicolinae because of the uncertain phylogenetic placement of putative early arvicoline fossils (Repenning, 1987;Fejfar et al., 2011).Thus we hypothesize an origin of arvicolines c. 6.4 Ma, with dentally distinct arvicoline rodents evolving slightly later.For ages of major tribes and clades of arvicolines found in our study, see Table 1.

CONCLUSIONS
A consolidated understanding of the phylogeny of arvicoline rodents has been warranted given their remarkable evolutionary history and abundance across high latitudes.We provide phylogenetic support for systematic hypotheses across Arvicolinae and some direction for future systematic and taxonomic work.We show that the earliest diverging arvicolines likely includes the Tribes Prometheomyini, Ondatrini, Lemmini, Dicrostonychini, and Pliophenacomyini.Prometheomys is probably the sister taxon of all other arvicolines, but work still needs to be done to solidify that hypothesis.The monophyly of bog lemmings (Synaptomys) is doubted, however their close relationship to the ''true lemmings'' (Myopus and Lemmus) is clear.Hyperacrius and Lemmiscus need to be sampled using phylogenomic data to better clarify their systematic positon.Clethrionomyini is a valid Tribe consisting of Caryomys, Eothenomys, Anteliomys, Craseomys, Alticola, and Clethrionomys.The paraphyly of Clethrionomys relative to Alticola needs to be explored.Pliomyini, Ellobiusini, Lagurini, Arvicolini, and Microtini represent a clade that forms a large portion of the arvicoline tree.We found some evidence for a monophyletic endemic North American Microtus, however this was not found in all analyses, and could benefit from larger scale sampling across the genome.Finally, we estimated divergence times among the major clades, which were concordant with the published fossil record.This will provide valuable insight into the evolutionary and paleobiogeographical history of this clade.Overall our understanding of the evolutionary history of arvicolines is improving and this study is a step in the right direction for both taxonomic and systematic clarity for Arvicolinae.
Lab at The University of Texas at Austin whose support was much appreciated.Having not generated any new molecular data for these analyses, we are immensely grateful to the many researchers whose work made this synthesis possible (see Appendix S1A for references to all authors who contributed data used in this study).We are grateful to Kenneth De Baets (editor), Connor J. Burgin (reviewer) and two anonymous reviewers for providing suggestions that greatly improved the paper.

Figure 1 Figure 2 Figure 3
Figure 1 Higher-order mitochondrial (Cytb and COI ) only analyses.(A) Maximum likelihood (ML) tree of mitochondrial dataset comprising 146 arvicolines and three members of Cricetidae.Tree was rooted with Neotoma fuscipes.For a species level tree see Appendix S1D.(B) Majority-rule consensus tree produced using Bayesian Inference (BI) methods of the mitochondrial only dataset comprising 146 arvicolines and three members of Cricetidae.For a species level tree see Appendix S1G.Full-size DOI: 10.7717/peerj.16693/fig-1

Figure 4
Figure 4 Time-calibrated phylogeny of Arvicolinae.Blue vertical dashed lines represent boundaries between geologic epoch.Red horizontal bar at each node represent the 95% HPD for the age of the node.Age in scale bar is in millions of years.PLE =Pleistocene.Orange Box =Holocene.For a species level timecalibrated phylogeny see Appendix S1J.Full-size DOI: 10.7717/peerj.16693/fig-4

Table 1 Divergence ages of various clades within this study.
Age is in millions of years before present.