Range contraction and expansion from glaciation has led to genetic divergence that may be particularly pronounced in fossorial species with low dispersal. The plains pocket gopher (Geomys bursarius) is a fossorial species that ranges widely across North America but has a poorly understood phylogeny. We used mitogenomes (14,996 base pairs) from 56 individuals across seven subspecies, plus two outgroup species, to assess genetic divergence from minimum spanning trees, measure genetic distances, and infer phylogenetic trees using BEAST. We found G. b. wisconsinensis was monophyletic with recent divergence. Further assessment is needed for G. b. major because it was paraphyletic and exhibited inconsistent groupings with other clades. Importantly, we identified G. b. illinoensis as being genetically distinct and monophyletic likely due to a unique colonization event eastward across the Mississippi River. Because G. b. illinoensis faces continued pressures from niche reduction and habitat loss, we recommend that G. b. illinoensis be considered an evolutionary significant unit warranting conservation actions to promote connectivity and restore suitable habitat. Such conservation efforts should benefit other grassland species including those originating from clades west of the Mississippi River that may also be evolutionary significant units.

We identified G. bursarius samples analyzed by Alexander et al. (2024) that had high DNA content based on microsatellite amplification to increase the likelihood of reconstructing complete mitogenomes. These samples included 41 toes from museum specimens from the Illinois Natural History Survey collected between 1934 and 1985, and 12 tissue samples from G. bursarius illinoensis live-trapped in 2018-2019, following appropriate guidelines (UIUC IACUC #17190). We sampled seven G. bursarius subspecies from across the range, excluding only G. b. ozarkensis. We included three samples from other Geomys species (one G. breviceps dutcheri, and two Geomys jugossicularis) as outgroups. Geomys jugossicularis was previously considered within the G. bursarius clade but was determined to be an independent species.

Arbor Biosciences (arborbiosci.com; Ann Arbor, MI, USA) completed library preparation, mtDNA-bait enrichment and capture, and genome sequencing. Arbor Sciences created customized mtDNA-baits using a myBaits® kit to enrich for Geomys mtDNA prior to sequencing. Baits were customized using a G. personatus reference mitogenome (TK24928, 16,817bp in length; Greenia, 2020). To account for potential DNA degradation in the 44 museum samples, we used 4-plex reactions with dual-round enrichment and subsequent sequencing of ~3 million paired-end reads. The 12 live-trapped individuals with non-degraded DNA were enriched in a single-round using a 12-plex reaction with subsequent sequencing of ~1 million paired-end reads. Enriched libraries were sequenced on the Illumina® NovaSeq platform to produce paired-end reads that were 150 base-pairs in length.

We used AdapterRemoval v.2.3.1 (Schubert et al. 2016) to trim adapters and consecutive Ns from the 5’ and 3’ ends of our sequencing reads. We used Burrow-Wheeler Alignment with maximal exact matches (BWA-MEM) v. 0.7.17 (Li and Durbin 2010) to align samples to the G. personatus mitogenome (G. personatus TK24928, 16,817bp in length; Greenia 2020). We used SAMTools v. 1.12 (Li et al. 2009) to remove duplicate and unmapped reads and reads with an alignment quality <30, then sorted and converted the alignment files to Binary Alignment Map (BAM) format. We calculated the breadth of coverage (i.e., the percentage of the genome that had ≥1 X-fold read coverage) and the average depth coverage (the average X-fold number of reads that mapped at any location across the genome).

We then imported the BAM files into Geneious Prime v. 2022.2.1 (Biomatters, Ltd www.geneious.com) and created a mitogenome consensus sequence that contained sites with a minimum depth of 5X reads using a majority base-call approach. In Mega v. 11.03 (Tamura et al. 2021), we manually curated our multiple-pairwise alignment by removing positions with high disagreement and low alignment quality that included portions of the D-loop and cyt-b regions. Because of the data-trimming and subsequent removal of cyt-b, we also individually aligned each read to the COX1 and cyt-b regions to a G. bursarius partial mitogenome (GenBank accession MZ030793.1; Greenia, 2020). We focused on these regions as COX1 had the lowest adaptive/synonymous replacement (~0), whereas cyt-b is commonly used in phylogenetic reconstruction and also has low adaptive/synonymous replacement in Geomys (Greenia 2020). We generated three datasets: (1) mitogenome-wide consensus sequences that were manually trimmed to remove disagreements, (2) consensus sequences based on reads that were independently aligned to the COX1 gene, and (3) consensus sequences based on reads independently aligned to the cyt-b region.

Literature Cited

Greenia, H. S. 2020. Comparative mitogenomics of the rodent suborder Castorimorpha. Thesis, Tarleton State University, Stephensville.

Li, H. et al. 2009. The Sequence Alignment/Map format and SAMtools. Bioinformatics 25:2078–2079.

Li, H., and R. Durbin. 2009. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics 25:1754–1760.

Schubert, M., S. Lindgreen, and L. Orlando. 2016. AdapterRemoval v2: Rapid adapter trimming, identification, and read merging. BMC Research Notes 9:88.

Tamura, K., G. Stecher, and S. Kumar. 2021. MEGA11: Molecular Evolutionary Genetics Analysis Version 11. Molecular Biology and Evolution 38:3022–3027.