The genus Microtus comprises 62 species of voles, distributed throughout North America, Europe, and Asia [1]. Microtus is believed to have undergone rapid speciation and diversification, with all speciation events occurring within the past 4 million years [2, 3]. It has been suggested that some nominal species, such as M. pennsylvanicus, contain cryptic diversity [4]. Microtus is an important model system across multiple biological disciplines, including studies of adaptation (e.g., [5]), infectious disease (e.g., [6]), parental care (e.g., [7]), and population dynamics (reviewed in [8]).

Currently, assembled genomes for four Microtus species have been deposited in GenBank: two European species (M. agrestis and M. arvalis) and two North American species (M. ochrogaster; [9], and M. oeconomus). The present study provides resources for two additional species: M. richardsoni (NCBI:txid111840) and M. montanus (NCBI:txid88450).

The North American water vole (M. richardsoni) occupies a large, disjunct distribution in the Pacific Northwest of North America, with habitat in the Cascades Mountains and the Rocky Mountains, spanning from southern Canada into central Utah. Four subspecies are currently recognized: M. r. arvicoloides in the Cascades Mountains, M. r. richardsoni in the Canadian Rocky Mountains, M. r. macropus in the central Rocky Mountains and Wyoming, and M. r. myllodontus in Utah. M. richardsoni is adapted to a semiaquatic lifestyle, relying on alpine and subalpine streams for creating burrows and escaping predators [10]. Like other semiaquatic mammals (e.g., otters), adaptations to this lifestyle have probably been driven by natural selection [11–13]. Water voles are among the largest species of Microtus, and their frequent movement makes runways of trampled vegetation along streams [10, 14]. Unlike most other vole species, M. richardsoni does not appear to experience regular population boom and bust cycles, although population size in the species may correlate with levels of precipitation [15]. Despite being listed in the International Union for Conservation of Nature (IUCN) Red List [16] as being of ‘Least Concern’, the species is listed as ‘critically imperiled’ in the Wyoming Natural Diversity Database owing to its specific habitat requirements, which can be substantially degraded by livestock grazing [17].

The montane vole (M. montanus) is partially sympatric with M. richardsoni, and can be found throughout most of the water vole’s range, with the exception of the Canadian Rockies. However, M. montanus can be found farther south and east, including areas of California, Nevada, Colorado, Arizona, and New Mexico [18]. The species has been divided into 15 subspecies, including M. m. canescens in the Cascades Mountains, M. m. nasus in the central Rocky Mountains, and M. m. amosus in northern Utah. Notably, M. montanus does not exhibit a break in its range in the Columbia Basin; probably because it is not restricted to riparian areas like M. richardsoni. The species is IUCN Red-Listed as being of ‘Least Concern’. However, M. m. arizonicus is ‘endangered’ according to the New Mexico State Game Commission Regulation [18], and M. m. ricularis is ‘of concern’, owing to its small range and declining population size [19].

The rapid radiation of Microtus voles has hindered systematic classification, leading to multiple taxonomic revisions and conflicting phylogenetic analyses [1, 20, 21]. Consequently, both species boundaries and relationships between species are difficult to infer.

Genomic resources within Microtus will help resolve these questions, and resources have steadily increased in recent years. Here, we present two nuclear and one mitochondrial genome assembly for M. richardsoni, as well as a structural and functional annotation for one of the M. richardsoni genomes. The subspecific classifications and disjunct range of the species means that M. richardsoni has been included in multiple phylogeography studies in the Pacific Northwest [22–24]. These studies were based solely on mitochondrial DNA, and the results of analyses that investigated species limits and demographic history were limited to inferences that can be derived from a single gene tree. These genomic resources will provide a rich source of data to address these knowledge gaps and aid in future studies of adaptation.

We also present single nuclear and mitochondrial genome assemblies for M. montanus. Genomic resources in M. montanus will provide a wealth of data to assess subspecies boundaries, quantify gene flow among subspecies, and aid in conservation efforts of threatened subspecies. Genome-level comparisons are also made between the new genome assemblies and other Microtus genome assemblies to examine differences in assembly quality and repeat content.

Frozen tissue from a single M. r. arvicoloides individual collected from the southern Cascades Mountain range (JMS_292; 44.016667N, −121.750000E; [22]) was sent to Hudson Alpha (Huntsville, AL, USA) for high molecular weight DNA extraction and 10× Chromium library preparation [25]. In the 10× method, each extracted DNA fragment receives a different barcode before the fragment is sheared for library preparation. After sequencing, these barcodes are used to connect sequencing reads for a more contiguous assembly. After sequencing with a single run on an Illumina HiSeqX, the resulting 150-bp paired-end reads were input into Supernova v. 2.2.1 for de novo genome assembly with –maxreads=all [26].

Additional tissue was obtained from a single M. r. macropus individual collected from the northern Rocky Mountains (JMG_88; 46.333333N, −114.633333E; [22]). DNA was extracted using a Qiagen DNeasy Blood and Tissue Kit, and the DNA was sent for library preparation and sequencing by Iridian Genomes, Inc (Bethesda, MD, USA). Then, 150-bp paired-end reads were sequenced on two runs of an Illumina HiSeqX. Genome assembly was performed using two different de Bruijn graph-based programs, SOAPdenovo v. 2.04 (SOAP, RRID:SCR_000689) and Discovar de novo v. 52488 (Discovar assembler, RRID:SCR_016755) [27, 28]. For SOAPdenovo, quality trimming was performed using fastQC v. 0.11.8 (FastQC, RRID:SCR_014583) and Trimmomatic v. 0.38 (Trimmomatic, RRID:SCR_011848) with settings ILLUMINACLIP: 2:30:10, LEADING:3, TRAILING:3, SLIDINGWINDOW:4:15, and MINLEN:36 [29, 30]. SOAPdenovo assemblies were performed with settings max_rd_len=150, avg_ins=300, reverse_seq=0, asm_flags=3, rd_len_cutoff=150, rank=1, pair_num_cutoff=3, and map_len=32. SOAPdenovo was run with k-mer values of 63, 89, 95, and 101 based on analysis of optimal k-mer values in kmerGenie v. 1.7051 [31]. Raw reads were used as input for de novo genome assembly with Discovar, as recommended in the program documentation.

To provide the most contiguous assembly for M. richardsoni, a hybrid assembly was performed using the ARCS+LINKS pipeline [32, 33]. The ARCS+LINKS pipeline uses barcoding information from the 10× Chromium reads to scaffold the contigs from a separate genome assembly. Barcoded reads from M. r. arvicoloides were mapped to the M. r. macropus Discovar assembly with bwa mem v. 0.7.17 [34] before converting the mapped reads to BAM format and sorting with SAMTools v. 1.8 (SAMTOOLS, RRID:SCR_002105) [35]. ARCS v. 1.0.5 and LINKS v. 1.8.6 were then run with settings –s 98 –c 5 –l 0 –z 500 –d 0 –r 0.05 –m 50-10000 –e 30000 and –d 4000 –k 20 –l 5 –t 2 –a 0.3 –o 0 –a 0.3 –z 500, respectively.

As part of a separate project, a single M. montanus individual from Utah (UMNH:Mamm:30891; 38.19381N, −111.5824E) was misidentified as M. richardsoni. DNA was extracted from the sample using a Qiagen DNeasy Blood and Tissue Kit before being sent to the University of California Davis Genome Center for library preparation and sequencing. Then, 150-bp paired-end sequences were collected with a single shared run on an Illumina NovaSeq. Species identity was confirmed using the Barcode of Life Database (BOLD [36]). Reads were checked and trimmed for quality with fastQC and Trimmomatic as above, before mapping reads to the mitochondrial cytochrome oxidase I (COI) sequence of M. r. macropus [37] using bwa mem. The resulting mapped reads were converted to BAM format, sorted, and indexed with SAMTools. PCR duplicates were identified and removed with Picard v. 2.3.0 (Picard, RRID:SCR_006525) [38]. Resulting reads were then piled with SAMTools mpileup using base and mapping quality scores of 30, consensus sequences were generated with bcftools v. 1.3.1 (SAMtools/BCFtools, RRID:SCR_005227) [39], and consensus sequences were converted to fastq format using vcfutils with a minimum depth filter of 5 and maximum depth filter of 10,000 [35]. The resulting sequence was input into BOLD. Owing to low sequencing coverage, de novo genome assembly was not appropriate for M. montanus. To provide a preliminary genome sequence, a reference-guided genome assembly was performed with RaGOO v. 1.1 [40]. Raw reads were input into Discovar to generate an initial genome assembly. Misassembly correction was performed with RaGOO, using reads trimmed with the same settings as the M. r. macropus reads, and RaGOO was then used to scaffold the Discovar contigs onto the M. r. arvicoloides assembly; this is more closely related to M. montanus than the other available Microtus genome assemblies [3]. Since M. montanus has less than half the chromosomes of M. richardsoni (2n = 22–24 in M. montanus versus 56 in M. richardsoni [41]), the possibility of structural errors in the M. montanus assembly was examined by calculating the percentage of reads that mapped back to the assembly using bwa mem and bamtools v. 2.2.2 [42].

The final assemblies were submitted to GenBank [43], where screening was performed to identify any contamination, and contaminated scaffolds were removed. All assemblies were evaluated with QUAST v. 5.0.1 (QUAST, RRID:SCR_001228) [44], bbmap v. 38.35 (BBmap, RRID:SCR_016965) [45], custom Python scripts [46], and BUSCO v. 5.0.0 (BUSCO, RRID:SCR_015008) using the Euarchontoglires reference set [47]. After comparing assembly statistics from the different assemblies of M. r. macropus, the Discovar assembly was selected as the best one because it had less fragmentation, higher N50 and lower L50 values, and a higher BUSCO score than the SOAPdenovo assemblies (Table 1). Genome sequencing of M. r. arvicoloides produced over 800 million reads and 47× genome sequencing coverage. The final genome assembly consisted of ∼32,000 scaffolds, with an N50 of 2.3 Mb (megabase pairs), 1.3% missing data (N), and a BUSCO score of 85.8%. Supernova estimated the length of the genome assembled to be ∼2.4 Gb (gigabase pairs) and the total genome size to be ∼2.6 Gb. M. r. macropus sequencing produced over 600 million reads and 35× coverage. Genome assembly with Discovar resulted in ∼1.6 million scaffolds, with an N50 of 16 Kb (kilobase pairs), 0.06% N, and a BUSCO score of 54.5%. Given the many programs available to perform de novo genome assembly from short reads, another program might have produced a more contiguous M. r. macropus assembly, but previous studies have shown that Discovar performs well compared to other programs [48, 49]. The hybrid assembly produced with the ARCS+LINKS pipeline had ∼1.6 million scaffolds, an N50 of 38 Kb, 0.09% N, and a BUSCO score of 59.8%. Because the hybrid assembly was of poor quality, it was not used for further analyses, and the M. richardsoni subspecies assemblies were kept separate. High fragmentation of the hybrid assembly is probably associated with fragmentation of the Discovar input assembly. Published results with this hybrid pipeline often include a much higher sequencing coverage of the input contigs to produce a better starting point for the pipeline. Therefore, in the future, additional Illumina sequencing with M. r. macropus might substantially improve the hybrid assembly. To produce the preliminary M. montanus genome, 108 million reads (13× coverage) were used, resulting in ∼13,000 scaffolds, an N50 of ∼3.1 Mb, 8.8% N, and a BUSCO score of 82.6%. Additionally, 89.3% of reads mapped back to the M. montanus assembly.

The complete mitochondrial genomes of M. r. arvicoloides and M. montanus were assembled using the genomic sequencing reads. Mitochondrial genomes were assembled by both mapping reads to a reference mitochondrial genome, and by using the reference-guided assembly program Novoplasty v. 4.1 [50]. For the mapping assembly, reads were mapped to the M. r. macropus mitochondrial genome using the same steps as in the M. montanus BOLD analysis. The mitochondrial assemblies were 16,285 bp and 16,268 bp in length, with an average depth of coverage of 7886× and 6805× for M. r. arvicoloides and M. montanus, respectively. Reference-guided mitochondrial assemblies with Novoplasty used the M. r. macropus mitochondrial genome as the reference, along with settings Genome Range  = 12,000–22,000, K-mer  = 33, Read Length  = 150, and Insert size  = 400. Because the M. r. arvicoloides dataset contained many reads, 25% of reads were subsampled for assembly, as suggested in the program documentation. The assemblies for M. r. arvicoloides and M. montanus were 16,298 bp and 16,319 bp in length, with average depths of coverage of 5131× and 14,713×, respectively.

To compare mitochondrial assemblies between methods, the assemblies were aligned using the MUSCLE plugin in Geneious v. R9 (Geneious, RRID:SCR_010519) with eight iterations and an open gap score of –1 [51, 52]. This comparison showed that the Novoplasty assemblies contained multiple insertions compared with the mapped assemblies and the reference mitochondrial genome. These insertions were up to 13-bp long in multiple genes, including trnT, trnK, and ATP8. Comparison with other Microtus mitochondrial genomes (M. ochrogaster; NC_027945.1 and M. fortis; NC_015243.1) showed that these insertions were only present in the Novoplasty mitochondrial assemblies. Therefore, the mapping assemblies were used for further analyses. The mapping assemblies for both species included ambiguous bases, which were much more frequent for M. montanus than M. r. arvicoloides. These could be the result of using the mitochondrial genome of a different subspecies (for M. r. arvicoloides) or species (for M. montanus) for mapping the reads. Additionally, the presence of nuclear DNA of mitochondrial origin (NUMTs [53, 54]) might have influenced these results. If mitochondrial segments have been incorporated into the nuclear genomes and subsequent mutations have occurred, both nuclear and mitochondrial sequences could be mapped to the same mitochondrial region during assembly, resulting in the ambiguous bases observed here. NUMTs are likely to be present, as they have been documented in other species of Microtus [55–57]. Both mitochondrial genomes were annotated using MITOS [58]. The annotations each consisted of 22 tRNA genes, two rRNA genes, and 13 protein-coding genes.

The available Microtus genome assemblies, M. agrestis (GCA_902806755.1), M. arvalis (GCA_007455615.1), M. ochrogaster (GCA_000317375.1), and M. oeconomus (GCA_007455595.1), were downloaded from GenBank. Assembly summary statistics were calculated using QUAST, bbmap, and custom Python scripts [46]. To compare repeat content among all genomes, including the three produced by the current study, repeats were first identified de novo using RepeatModeler v. 2.0.1 (RepeatModeler, RRID:SCR_015027) [59]. RepeatMasker v. 4.0.9 (RepeatMasker, RRID:SCR_012954) was then used to further identify repeats using a combined repeat library that included the repeats identified from RepeatModeler and those from the RepeatMasker Rodentia database [60]. The percentage of the genome comprising each type of repeat element was extracted from the RepeatMasker log file for each genome assembly.

All genome assemblies used some form of Illumina sequencing (Table 2), although assembly continuity varied greatly among assemblies, from 1366 scaffolds in M. agrestis to 1.6 million scaffolds in M. r. macropus. Genome coverage was similarly varied, from 13× in M. montanus, to 35× in M. r. macropus, to 77× in M. arvalis and M. oeconomus. The percentage of repetitive regions ranged from 31.7% in M. montanus to 44.1% in M. arvalis (Figure 1). Repeat content did not appear to be associated with phylogenetic relatedness, since repeats between the two subspecies of M. richardsoni were not more similar to each other than to other Microtus species. However, it is possible that the repeat content is affected by the continuity of the genome assemblies. Further research is needed to confirm this relationship.

The M. r. arvicoloides genome assembly was annotated with the MAKER v. 2.31.10 (MAKER, RRID:SCR_005309) pipeline [61], loosely following the tutorial by Daren Card [62]. Briefly, the pipeline comprises masking repeats followed by multiple rounds of annotation, with both evidence-based and ab initio gene models. Repeats were identified as described above. Complex repeats were then extracted from RepeatMasker results using grep with keywords “Satellite” and “rich”. Within Maker, the model_org argument was set to “simple” so Maker would soft-mask simple repeats, and the RepeatMasker results were provided to hard-mask complex repeats. Evidence-based gene discovery used protein and mRNA sequences from the previous genome annotation of M. ochrogaster (GCF_000317375.1), as well as an additional RNA-Seq assembly from M. pennsylvanicus (GSM3499528 [63]). Hidden Markov models (HMMs) for ab initio gene prediction were trained using both SNAP and Augustus v. 2.5.5 (Augustus: Gene Prediction, RRID:SCR_008417) [64, 65]. With SNAP, gene models identified by MAKER were filtered using an annotation edit distance (AED) of 0.5 and an amino acid length of 50. After validating these models with SNAP’s Fathom utility, removing likely errors, and including 1000 bp surrounding each training sequence, the training sequences were passed to the hmm-assembler script. For Augustus, training sequences plus 1000 bp on each side were obtained from the first round of MAKER mRNA annotations. Augustus was used to train the HMM using the –long option in BUSCO and the Euarchontoglires reference set. MAKER was then run again with the previously annotated gene models, and the HMM models from SNAP and Augustus. After the initial MAKER run, two cycles of ab initio gene prediction and annotation with MAKER were performed. To prevent overfitting, results were compared after each round of MAKER. Because the increase in AED score was minimal between the first and second rounds of ab initio gene prediction, further analysis was conducted on the results after the first round only. This round annotated ∼24,500 genes, with a mean gene length of 7445 bp (Table 3), which is within the range found in previous studies of M. ochrogaster (22,427 genes; GCF_000317375.1) and Arvicola amphibious (25,136 genes; GCF_903992535.1). Of these annotations, all occurred on scaffolds greater than 1 Kb in length, and 97% occurred on scaffolds greater than 10 Kb in length.

Functional annotation of the M. r. arvicoloides genome was performed using GOfeat [66], an online functional annotation tool that uses multiple protein databases including UniProt [67], InterPro[68], and Pfam [69]. An input file for GOfeat was generated by supplying the genome assembly FASTA file and the MAKER General Feature Format (GFF3) file to the Python package gffread v. 0.12.1 [70]. GOfeat annotated 83.49% of genes. Biological processes accounted for 42.46% of annotations, cellular components accounted for 30.29%, and molecular functions comprised 27.25%. The most frequent gene ontology (GO) terms were ‘positive regulation of transcription by RNA polymerase II’, ‘negative regulation of transcription by RNA polymerase II’, and ‘DNA-templated regulation of transcription’ for biological processes; ‘cytoplasm’ and ‘plasma membrane’ for cellular components, and ‘metal ion binding’ and ‘calcium ion binding’ for molecular functions.