Newly assembled pocket gopher genomes can facilitate conservation management of biodiversity

Abstract

Texas exhibits one of the richest levels of pocket gopher diversity in the United States. Three genera (Cratogeomys, Geomys, and Thomomys) and 11 species are found in Texas. It is not surprising given the diversity of the Texas landscape (ecoregions, life zones, substrates, and vegetation) that these species are further subdivided into 29 subspecies in Texas alone. Pocket gopher distributions are determined by availability of suitable soil types and therefore often occur in small, isolated populations. For some taxa, limited distribution and ultimately small deme sizes result in populations that may require attention from a regulatory and management perspective. For many Texas pocket gopher subspecies, insufficient information exists to make sound recommendations relative to conservation status and needs despite decades of research collecting and evaluating data based on morphometrics, distributions and habitat preferences, karyotypes, allozymes, and mitochondrial DNA. As such, there is precedent for elevating pocket gopher subspecies to species after evaluation of available data, as well as subsuming subspecies into a broader taxonomic group. We used genomic techniques to identify genetically defined operational taxonomic units (OTUs) of pocket gophers to improve knowledge and understanding of pocket gopher distributions within the state. Using tens of thousands of single nucleotide polymorphisms, we determined the number of OTUs in each genus to be 5 for Thomomys bottae subspecies, 8 for Geomys species, and 5 for Cratogeomys castanops subspecies in Texas. In general, these data agree with current taxonomic hypotheses regarding Geomys and C. castanops; however, many T. bottae groups present similar genetic patterns that do not merit subspecies status based on these data, suggesting a more conservative classification of T. bottae in Texas and southeastern New Mexico that could facilitate conservation efforts, should they be necessary.

Resumen

Texas presenta uno de los niveles de diversidad de tuzas (gofers de bolsillo) más altos de los Estados Unidos de América. En Texas se encuentran tres géneros, Cratogeomys, Geomys y Thomomys, y once especies. Esto no es sorprendente, dada la diversidad del paisaje de Texas (ecorregiones, zonas de vida, sustratos y vegetación), que estas especies estén subdivididas en 29 subespecies. Las distribuciones de las tuzas están determinadas por la disponibilidad de tipos adecuados de suelo y, por lo tanto, a menudo corrresponden a poblaciones pequeñas y aisladas. Para algunos taxones, la distribución limitada y, en última instancia, los grupos de pequeño tamaño resultan en poblaciones que pueden requerir atención desde una perspectiva regulatoria y de gestión. Para muchas subespecies de tuzas de Texas existe información insuficiente para hacer recomendaciones sólidas sobre el estado y las necesidades de conservación, a pesar de décadas de investigación que recopilan y evalúan datos basados en morfometría, distribuciones y preferencias de hábitat, cariotipos, alozimas y ADN mitocondrial. Así pues, existe un precedente para elevar las subespecies de tuzas a especies después de evaluar los datos disponibles, así como para incluir las subespecies en un grupo taxonómico más amplio. Utilizamos técnicas genómicas para identificar unidades operativas de taxonomía definidas genéticamente (OTUs) de tuzas para mejorar el conocimiento y la comprensión de las distribuciones de dentro del estado. Usando decenas de miles de polimorfismos de nucleótido único (SNPs), determinamos que el número de OTUs en cada género es 5, 8 y 5, para las subespecies de Thomomys bottae, las especies de Geomys y las subespecies de Cratogeomys castanops, respectivamente, en Texas. En general, estos datos coinciden con las hipótesis taxonómicas actuales sobre Geomys y C. castanops; sin embargo, muchos grupos de T. bottae presentan patrones genéticos similares que no justifican el estatus de subespecie basado en estos datos, lo que sugiere una clasificación más conservadora de T. bottae en Texas y el sureste de Nuevo México que podría facilitar los esfuerzos de conservación, si fueran necesarios.

A species is the unit of currency in biodiversity conservation (Wilson 2017; Coates et al. 2018), yet defining a “species” often is not a simple task. Conservation managers, under the direction and influence of societal, legal, and scientific stakeholders, are tasked with determining which groups of organisms are deserving of special protection. But, as mentioned above, “species” means different things to different governing bodies and many of them recognize infraspecific designations for possible listing candidacy. As the name suggests, subspecies are a classification below the level of species that arose in the 19th century to describe variation within a species (Wilson and Brown 1953; Durrant 1955; Lidicker 1962). As such, subspecies can and are defined by any manner of criteria, including morphology, ecology, geography, behavior, and/or genetics. Because speciation is rarely instantaneous (Hendry et al. 2007), older named subspecies may eventually represent incipient species (Mayr 1982). Members of infraspecific groups, if threatened with small population size, habitat destruction, or demographic stochasticity, can therefore be listed and offered protections. However, the opposite may be true—perhaps an infraspecific group, thought to exhibit “evolutionary significance,” may receive unneeded protections (Frankham et al. 2012). Defining or delimiting subspecies is also challenging, and in the genomics era, this rank has come under scrutiny (see Burbrink et al. 2022).

A group whose taxonomy has changed considerably over the years is that of the well-studied North American pocket gophers (family Geomyidae). Pocket gophers are subterranean rodents with interesting ecological, physiological, and genetic characteristics. Vast morphological diversity is displayed across the 6 genera of pocket gophers, and as such there are hundreds of described subspecies (Nowak and Walker 1999). Morphological variation in body size, pelage coloration, and cranial characters seen in this group of rodents has been hypothesized to have different causes. For example, groups that, on average, reach a larger body size may have better nutrition (e.g., alfalfa fields; Patton and Brylski 1987); pelage coloration is similar to the excavated soils in which they are found, suggestive of camouflage (Ingles 1950; Krupa and Geluso 2000); or that incisor procumbence is a function of soil hardness (Lane 1965; Marcy et al. 2016). With such variation in phenotype existing in pocket gophers, it is challenging to ascertain if these differences are changes due to direct local environmental pressures, genetic drift, or some other biological factors. Because of the high phenotypic variation in morphology, it becomes important to analyze and group these rodents into units based on genetic cohesiveness rather than on morphological characteristics alone.

Although many subspecies are found in California, Arizona, and New Mexico, a pocket gopher “hotspot” is the state of Texas (Nowak and Walker 1999). Pocket gophers are distributed throughout the state but tend to occur in highly fragmented populations due to varying habitat requirements, namely soil suitability (Miller 1964). Because of this patchy distribution and complexity of convergent evolution and in some cases sampling difficulties, the status of 32 taxa comprising the genera Geomys, Thomomys, and Cratogeomys that occur in Texas remains ever-changing. It is possible that some of these small and isolated groups may be at an elevated conservation risk (Nevo 1979). As such, the main research aim of this paper is to use genomic methods, including assembly of 2 new gopher genomes, to delimit operational taxonomic units (OTUs) that provide an alternative and perhaps more digestible view on Texas and southwestern New Mexico pocket gopher taxonomy in hope of maintaining biodiversity of this group. We also highlight pocket gopher groups that we could consider prioritizing for additional surveys in the near future.

Materials and methods

Fieldwork

We surveyed many accessible locations in and around the type localities listed in Table 1 and shown in Fig. 1 and caught gophers in these areas (Appendix  I). Once an animal was captured, a series of standard measurements were taken of the animal, and then the animal was euthanized using isofluorane. Extracted tissues (heart, kidney, muscle, liver, lung, and spleen) were prepared for deposition in Natural Science Research Laboratory (NSRL) at the Museum of Texas Tech in Lubbock, Texas. Individuals were collected in adherence to the Texas Tech University’s Institutional Animal Care and Use Committee guidelines (#20002-01), and those guidelines set forth for research involving wild mammals found in Sikes et al. (2016).

Genome sequencing and assembly

DNA extraction and long-read sequencing were performed at the Brigham Young University DNA Sequencing Center. DNA was extracted from liver subsamples from G. bursarius (NSRL TK 200958) and T. bottae (NSRL TK197808) using the Qiagen Genomic-tip 10/G kit (Qiagen, Germantown, Maryland). Sequencing was performed using a PacBio Sequel II on 2 SMRT cells per individual. Libraries were prepared using the SMRTbell Express Template Prep Kit 2.0 (Pacific Biosciences, Menlo Park, California) according to the manufacturer’s protocol.

Short read sequencing was performed at Novogene, Inc. DNA was extracted using the Qiagen DNeasey prep kit. Genomic DNA was randomly fragmented into 350-bp, end-repaired, adenylated, ligated with Illumina sequencing adapters, and PCR-enriched. Final libraries were purified (AMPure XP system) followed by quality control, including size verification by Agilent 2100 Bioanalyzer (Agilent Technologies, California) and molar concentration using real-time PCR.

Preassembly estimates of genome sizes were determined using k-mer counts and Jellyfish (Marcais and Kingsford 2011) and GenomeScope (Vurture et al. 2017) with k-mer values ranging from 25 to 65. Initial genome assembly was accomplished with Flye (Kolmogorov et al. 2019) and followed by polishing with Illumina reads using POLCA (Zimin and Salzberg 2020). Final assemblies were evaluated with BUSCO v4.1.4 (Simao et al. 2015; Seppey et al. 2019) using the vertebrate and mammalian protein sets (vertebrate_odb10 and mammalia_odb10).

DNA extraction and genomic sequencing

For all samples, DNA was extracted using the Qiagen DNeasy Blood and Tissue spin column protocol (Qiagen, Venlo, The Netherlands) following the manufacturer’s protocol. However, instead of an elution volume of 200 µL, DNA was eluted to final volumes between 100 and 125 µL to increase sample concentration. DNA concentration was fluorometrically quantified using the Qubit 3.0 broad-range assay (Invitrogen, Life Technologies, Carlsbad, California).

Prelibrary samples were sent to Admera Health in South Plainfield, New Jersey, United States, for double-digest restriction site-associated DNA sequencing or ddRADseq (Baird 2008; Peterson et al. 2012). All dual-indexed libraries were pooled into a single final library that was sequenced across multiple Illumina HiSeq 2500 lanes.

After receiving Illumina data from the sequencing center, subpar reads were filtered out of the data set using the AfterQC “after.py” pipeline (Chen et al. 2017). Subpar reads were characterized as having low-quality scores (PHRED score < 15), mismatched reads, short reads (<35 base pairs), too many ambiguous nucleotides (greater than 40% of read), or excessive homopolymer regions. If a read satisfied one of the above criteria, it was removed from the data. Retained reads were aligned to a draft T. bottae genome or draft G. bursarius genome using the “bwa-mem” module within the Burrows–Wheeler aligner (Li and Durbin 2009). Cratogeomys was mapped to the Geomys genome. This procedure improved mapping quality of our generated RAD-Seq reads.

SNP-based population and phylogenetic analyses

Genotypes were analyzed and exported to GenePOP and STRUCTURE using Stacks v2.01 (Catchen et al. 2013). The “populations” module in Stacks was used iteratively to identify the optimized data set whereas adjusting certain parameters such as minimum number of populations (-p) and minimum samples per population (-r) flags (Catchen et al. 2013; Mastretta-Yanes et al. 2015). Only one random single nucleotide polymorphism (SNP) per locus was analyzed to minimize possible confounding effects of covariation between SNPs.

STRUCTURE, which incorporates a Bayesian algorithm (Pritchard et al. 2000), was used to assign individual likelihood of belonging to a population. A random subset of at least 10% of the final data was used to increase processing speed. To test model congruence, SNPs were randomly selected 5 times with 25,000 burn-in iterations and an additional 75,000 Markov Chain Monte Carlo repetitions in STRUCTURE at each k, which ranged from 1 to 15. Output from STRUCTURE was curated using STRUCTURE HARVESTER (Earl and vonHoldt 2012) and CLUMPP (Jakobsson and Rosenberg 2007). STRUCTURE HARVESTER was also used to find the optimal value of k for the data set using the Evanno method (Evanno et al. 2005). We used Distruct for final visualization of STRUCTURE plots (Rosenberg 2004).

Additionally, Bayesian phylogenetic methods were used to infer relationships between genera using a concatenation of SNPs. The Bayesian tree was inferred using MrBayes, incorporating the gamma rate and the Kimura 2 parameter, which is a simple model that is widely used with closely related sequences (Ronquist et al. 2012).

Determination of OTUs

Genetically determined OTUs were determined using the following criteria: (i) an F_ST value equal to or greater than 0.25 (Takahata and Nei 1984; Meirmans and Hedrick 2011); (ii) after the most likely value of k for the group was selected, the proportion of membership for a particular cluster was greater than 0.5 using STRUCTURE (simple majority); and (iii) in Bayesian phylogenetic reconstruction, clades were formed based on 99% similarity (1% difference). These criteria are borrowed from studies of microbial diversity, where OTUs are grouped based on the de facto 97% identity (or 3% difference).

The average number of clades determined from all analyses was used for the final number of clades per genus. Each genus was evaluated on a case-by-case basis.

Results

Genome assembly statistics

PacBio reads numbered 15.3 M and 14.1 M, and Illumina reads numbered 421.6 M and 537.1 M (PE 150) for G. bursarius and T. bottae, respectively. PacBio reads averaged 21,174 bp and 22,035 bp for G. bursarius and T. bottae, respectively. Preassembly estimates of genome size ranged from 1.84 to 2.08 for G. bursarius and 1.98 to 2.51 for T. bottae and were very similar to final assembly sizes (Table 2). The BUSCO analysis indicated excellent representation of mammalian genes, with 92.4% of complete mammalian BUSCOs identified in both genomes and low rates of fragmented and missing orthologs (Supplementary Data SD1). Both assemblies have been deposited in the NCBI database (accession numbers SAMN31140755 for Thomomys and SAMN31140756 for Geomys).

Thomomys bottae

We analyzed 14 named subspecies of T. bottae from southeastern New Mexico and southwestern Texas (Table 1). An average of 158,377 loci were generated for each subspecies, and a total of 35,909 SNPs were evaluated from these (Supplementary Data SD2).

F_ST values were generally under 0.25 for pairwise comparisons (Supplementary Data SD3). The 2 subspecies with the greatest pairwise F_ST were T. b. lachuguilla and T. b. tularosae at 0.504. According to the STRUCTURE HARVESTER, the most likely value of k for the Thomomys data set was 5. Except for the New Mexico subspecies, T. b. tularosae, all groups show admixture (Fig. 2). Whereas most groups show membership coefficients to at least 1 cluster greater than 0.5, T. b. confinalis, T. b. limitarus, and T. b. robertbakeri show almost equal membership coefficients in 2 clusters. It is clear from the Bayesian phylogeny (Fig. 2) that New Mexico subspecies T. b. ruidosae and T. b. tularosae are divergent from all T. bottae subspecies examined. Remaining T. bottae subspecies in Texas show little genetic differentiation with the exception of T. b. actuosus.

Geomys species and subspecies

We examined 18 groups in southwest Texas (Table 1). An average of 161,907 loci were generated for each subspecies, and a total of 79,651 SNPs were evaluated from these (Supplementary Data SD4).

F_ST values were larger in this genus than in Thomomys and Cratogeomys, likely because comparisons were often between recognized species, not subspecies. The greatest pairwise F_ST value was between G. arenarius brevirostris and G. tropicalis at 0.811 (Supplementary Data SD5). STRUCTURE and STRUCTURE HARVESTER determined the most likely value of k to be 7 (Fig. 3). Groups with moderate levels of admixture were G. b. major and G. streckeri. It is evident from the STRUCTURE plots that some subspecies have been assigned 1 taxonomic group a priori, but the taxonomy was not current (Baker and Glass 1951). Two such cases are G. b. dutcheri and G. b. brazensis, which show similar memberships to that of G. breviceps, and G. b. ammophilus that shows similar memberships to G. attwateri.

The Bayesian tree (Fig. 3) constructed for Geomys groups examined reveal 2 clades: species G. arenarius, G. bursarius, G. juggosicularis, G. knoxjonesi, and G. texensis sharing common ancestry. Geomys attwateri, G. breviceps, G. personatus, G. streckeri, and G. tropicalis form the second clade.

Cratogeomys castanops

Seven groups of Cratogeomys castenops were examined in southwest Texas and southeastern New Mexico (Table 1). An average of 144,607 loci were generated for each subspecies and from these, a total of 33,284 SNPs were examined (Supplementary Data SD6). For pairwise F_ST comparisons, both C. c. angusticeps and C. c. tamaulipensis show moderate levels of genetic differentiation, with the greatest pairwise comparison of 0.501 identified by comparing C. c. tamaulipensis to C. c. perplanus (Supplementary Data SD7). Only C. c. clarkii had all values below 0.25, indicating genetic similarity between that subspecies and all other groups examined. The most likely value of k via STRUCTURE and STRUCTURE HARVESTER was 5. Cratogeomys c. angusticeps, C. c. parviceps, and C. c. perplanus show no admixture with other groups, whereas C. c. clarkii, C. c. dalquesti, C. c. tamaulipensis, and C. c. lacrimalis do show evidence of admixture (Fig. 4). The Bayesian phylogeny (Fig. 4) supports C. c. angusticeps and C. c. tamaulipensis as sister taxa to all other C. castanops evaluated. The tree also supports C. c. clarkii, C. c. parviceps, C. c. lacrimalis, and C. c. dalquesti sharing common ancestry, with C. c. perplanus as their sister taxon.

Discussion

We believe that it is fortunate that a large-scale analysis of Texas pocket gophers was conducted at this time. Many groups have never been evaluated at genomic-scale resolution. Results herein add to the current taxonomy and could influence management efforts, especially groups highlighted in Table 3. Research into many of these same groups using the mitochondrial Cytochrome b gene to update taxonomic recommendations was recently published (Bradley et al. 2023).

Thomomys bottae (=Thomomys baileyi; Bradley et al. 2023)

Based on population and phylogenetic analyses, we suggest 5 clades for T. bottae in Texas and southeastern New Mexico (Supplementary Data SD8). Bailey’s assessment of Thomomys used morphology as the key determinant of species (Bailey 1915). As both evolutionary theory and scientific instruments have changed drastically in the last 100 years, so has the T. bottae–umbrinus species complex (Álvarez-Castañeda 2010). Within the last 2 decades, there have been at least 2 major genetic assessments of T. bottae subspecies from the region under study that have led to vastly different interpretations of taxonomic classifications in the group. Wickliffe (2004) concluded that although there was great morphological variation, these groups displayed little genetic divergence in the 7 taxa studied. On the other hand, Álvarez-Castañeda (2010) examined the Cytochrome b gene and recommended dividing the T. bottae–umbrinus complex into 8 phylogenetically distinct species. Obvious reasons for disagreement between these 2 genetic investigations could be related to the data set used, species concept applied, or underlying biological factors, such as cytonuclear discordance. For these reasons, this independent study of pocket gophers in this region using many noncoding loci is particularly appropriate.

Thomomys bottae ruidosae–tularosae group

All lines of evidence point to T. b. ruidosae and T. b. tularosae being diverged from the rest of T. bottae, though T. b. ruidosae may harbor alleles from other T. bottae subspecies in Texas (Fig. 2). It is possible that T. b. tularosae is a distinct species, and T. b. ruidosae is a hybrid form. Samples sizes were low for T. b. ruidosae in this study, so further investigation into this speculation is warranted. T. umbrinus, another pocket gopher in subgenus Megascaphus was shown by Patton (1973) to hybridize with T. bottae in the Patagonia Mountains of Arizona, so this supposition is not entirely out of the question. Unfortunately, our sampling lacked any individuals outside of T. bottae in southwestern Texas and southeastern New Mexico. Nevertheless, a question that could be easily followed up and answered is the placement of T. b. tularosae and T. b. ruidosae as either T. bottae (=baileyi), as it is currently, or within T. umbrinus. Keeping these 2 subspecies with T. bottae is supported by Beauchamp-Martin et al. (2019), though Bradley et al. (2023) considers elevating T. b. ruidosae to species status based on genetic divergence of the Cytochrome b gene.

Thomomys bottae confinalis–limitarus–robertbakeri group

Bayesian trees support monophyly of T. b. limitarus, T. b. robertbakeri, and T. b. confinalis. These groups occur in the southern Trans-Pecos and extend east to the Edwards Plateau. The easternmost subspecies of T. bottae in Texas is T. b. confinalis, which can be found as far east as Mason County, Texas. In 2019, a new taxon was suggested by Beauchamp-Martin et al. (2019) based on cranial morphological characters. The authors argued that there is no size-related natural diagnostic cutoff in the region occupied by T. b. limitarus and T. b. confinalis, suggesting that the new taxon would solve this. However, other mechanisms such as phenotypic plasticity and genetic drift could be responsible for the distribution in size in this group. In the phylogenetic trees inferred here, T. b. robertbakeri (the new taxon) and T. b. confinalis (the taxon whose demarcation was reduced in Beauchamp-Martin et al. 2019) show an unresolved relationship (Fig. 2). STRUCTURE analysis demonstrates that T. b. confinalis, T. b. limitarus, and T. b. robertbakeri all show similar assignments, and therefore are placed in the T. b. confinalis–limitarus–robertbakeri group. This assignment of retaining T. b. confinalis and T. b. robertbakeri based on size (Beauchamp-Martin et al. 2019) and subsuming T. b. limitarus would be consistent with our data. Thomomys b. limitarus may also belong to the T. b. pervarius–spatiosus or the T. b. confinalis–limitarus–robertbakeri clade, but further investigation is needed to delimit the boundaries of these populations. Such investigation could be hindered by the encroachment of C. castanops into the Trans-Pecos (Reichman and Baker 1972).

Thomomys bottae limpiae–texensis group

Thomomys bottae limpiae and T. b. texensis make up the T. b. limpiae–texensis group. The type localities for these populations are in the same county (Jeff Davis County in Texas), so it is unsurprising that clustering analysis would support grouping them into a single clade. According to STRUCTURE, T. b. limpiae shows greater probability of assignment with the lowland groups (T. b. pervarius–spatiosus and T. b. confinalis–limitarus–robertbakeri) than does T. b. texensis (Fig. 2), a result that is supported by the type localities of each group (Table 1). Combining T. b. texensis and T. b. limpiae based on size (Beauchamp-Martin et al. 2019) would be consistent with the data presented here.

Thomomys bottae pervarius–spatiosus group

Perhaps the most elusive group is the T. b. pervarius–spatiosus populations because efforts to capture Thomomys in these areas were rather unsuccessful. Type localities for both groups are in the lowlands in the Trans-Pecos, Shafter, Texas; and Alpine, Texas, respectively. Beauchamp-Martin et al. (2019) suggested aligning T. b. pervarius with T. b. lachuguilla and that T. b. spatiosus is unique. Lack of capture success for T. b. spatiosus and T. b. pervarius limited the sample size for this group, and it is possible that these individuals could be included in another group, as demonstrated by the phylogenetic analysis (Fig. 2).

Thomomys bottae guadalupensis–lachuguilla–actuosus–scotophilus–pectoralis group

All Thomomys populations found on mountain ranges on the Franklin Mountains and Guadalupe Mountains, which are north of the Davis Mountains, form this clade. Historically, many of these groups received subspecies designation because of geography (Jones and Baxter 2004), but STRUCTURE analysis shows little to no genetic differentiation (Fig. 2). This result suggests high levels of gene flow among all groups, which would have to extend into the lowlands between mountain ranges. Additionally, soil texture is possibly too rocky to support the larger-bodied Cratogeomys and may not be suitable for human settlement, in which case may influence collection effort and accessibility (i.e., roads). Phylogenetic analysis suggests that T. b. actuosus is most divergent from the other 4 and Beauchamp-Martin et al. (2019) classified T. b. actuosus as its own subspecies.

Similar to Beauchamp-Martin et al. (2019), we suggest synonymizing some T. bottae groups in Texas. However, specific groupings between Beauchamp-Martin et al. (2019) and this study are discordant, and it is not uncommon for morphology and different types of molecular data to disagree (Wickliffe et al. 2004; Álvarez-Castañeda et al. 2010), especially when comparing data across the nuclear genome to 1 mitochondrial locus or morphological character. For example, T. b. actuosus forms its own clade in Beauchamp-Martin et al. (2019), whereas in this analysis, T. b. actuosus forms this group clade with 4 other subspecies—T. b. lachuguilla, T. b. scotophilus, T. b. pectoralis, and T. b. guadalupensis. The latter 3 in Beauchamp-Martin et al. (2019) are grouped together with T. b. texensis, whereas in the present study T. b. texensis groups with T. b. limpiae. Genetic results herein indicate that T. b. lachuguilla is similar to others in this group even though it is found along the Rio Grande corridor. Bradley et al. (2023) suggested subsuming T. b. lachuguilla with T. b. limitarus (in part) and T. b. pervarius (in part), which is different from what Beauchamp-Martin et al. (2019) suggested, by placing T. b. limitarus (in part) with T. b. lachuguilla and (in part) with T. b. robertbakeri. These differences indicate a group in need of further study. Finally, results here place T. b. limitarus in the T. b. confinalis–limitarus–robertbakeri clade. Thomomys b. baileyi was not included in this study because of the lack of genetic material.

Geomys spp

Studies by Block and Zimmerman (1991), Jolley et al. (2000), Sudman et al. (2006), Genoways et al. (2008), and Chambers et al. (2009) evaluated many of the pocket gophers presented here using a variety of molecular markers. In this study, at least 3 subspecies, G. b. ammophilus, G. b. brazensis, and G. juggosicularis (Baker and Glass 1951; Honeycutt and Schmidly 1979; Bohlin and Zimmerman 1982) classified as their own taxon have genetic composition and distributions similar to another defined taxon. Results generally concur with Chambers et al. (2009) which used the nuclear-encoded interphotoreceptor retinoid-binding protein gene (Rbp3) and mitochondrial 12S ribosomal RNA (12S rRNA) genes to delimit 4 species groups: G. bursarius, G. breviceps, G. personatus, and G. pinetis—the latter, which is found in the southeastern United States (Pembleton and Williams 1978), was not included in the present work.

North and West Texas groups: Geomys arenarius, G. bursarius, G. knoxjonesi, G. texensis.

Geomys arenarius arenarius and G. a. brevirostris show similar genetic patterns and are clearly distinguished from other such groups. Therefore, these remain valid taxa. Geomys bursarius is a widespread species, with populations extending into Canada. Within Texas, this species can be found in the north-central part of the state, the high plains, and into the extreme northwest panhandle (Supplementary Data SD9). Jolley et al. (2000), Sudman et al. (2006), Genoways et al. (2008), and Chambers et al. (2009) used Cytochrome b data to argue that G. juggosicularis was divergent from G. bursarius and was aligned with gophers collected from the upper Midwest. In this study, G. b. major and G. juggosicularis show similar genetic assignments and affiliations (Fig. 3). The population found in the northern panhandle of Texas and identified as G. juggosicularis may belong be G. b. major. These 2 groups, then, form the G. b. major clade, echoing Villa and Hall’s (1947) placement of G. juggosicularis.

Geomys knoxjonesi, though it aligns with G. arenarius in the STRUCTURE data (Fig. 3), shows enough genetic dissimilarity from G. arenarius to retain its species status and confirms the placement by Sudman et al. (2006). Central Texas pocket gophers (G. texensis texensis, G. t. llanensis, and G. t. bakeri) are genetically similar to each other yet distinct from all others in all analyses and so comprise valid taxa, corroborating the results of Sudman et al. (2006) and Chambers et al. (2009).

South and East Texas groups: Geomys attwateri, G. breviceps, G. personatus, G. streckeri.

Geomys attwateri and G. b. ammophilus align with each other in STRUCTURE and phylogenetic analyses (Fig. 3), and accordingly, Bradley et al. (2023) recommend subsuming G. b. ammophilus with G. attwateri. These pocket gophers are genetically distinct from G. breviceps sagittalis, G. b. dutcheri, and G. b. brazensis, which fall under G. breviceps. In Bradley et al. (2023), G. b. brazensis had been suggested to be named G. brazensis brazensis to distinguish it from G. breviceps. Furthermore, the new recommendation is to rename Geomys b. dutcheri to Geomys breviceps dutcheri. Finally, interpretation from the results from Bradley et al. (2023) suggests that G. breviceps sagittalis could, in part, be a new species. More work is needed here to assess the strength of that assessment.

The G. personatus groups have a complicated taxonomy, as outlined in Chambers et al. (2009). Based on our results, G. p. maritimus, G. p. davisi, G. p. megapotamus, and G. tropicalis all fall within the same clade in Texas. However, G. tropicalis is the most divergent member of Geomys based on chromosomes and should be treated as a distinct and valid taxon (Kim 1972). When compared to other populations within this group, G. p. personatus had no private alleles, which suggests that this group has likely been subsumed into one (or more) of the other groups or that the tissue used for genetic analysis (toe clips from museum samples) did not provide the same quality of DNA as liver samples. Sudman et al. (2006) noted that further study was needed to determine the status of G. b. jugossicularis, G. p. davisi, G. p. maritimus, and G. b. sagittalis. From this study, none of these groups appear genetically distinct enough to rise to the level of a species. The phylogenetic tree generated from these data also disagrees with both Sudman et al. (2006) and Chambers et al. (2009) with the assertion that if G. tropicalis is recognized as distinct from G. personatus that the latter must be paraphyletic—in Fig. 3, G. personatus is monophyletic with respect to G. tropicalis. Two subspecies were not captured during this study: G. p. fallax and G. p. fuscus. These were also not included in Sudman et al. (2006), a sobering indication that these groups may have become extinct. Their inclusion may have proved significant, possibly helping to resolve relationships within G. personatus. Finally, found in southern Texas counties of Dimmitt, Zavala, and Jim Hogg (Supplementary Data SD9), G. streckeri represents a valid taxon, with curious STRUCTURE results illustrating moderate levels of assignment with G. personatus, G. breviceps, and G. attwateri (Fig. 3).

Cratogeomys castanops

Cratogeomys castanops is a widespread member of the genus and was re-elevated from Pappogeomys in 1982 based on synapomorphic alleles using electrophoretic data from allozymes (Honeycutt and Williams 1982). Morphological synapomorphies are also apparent: a sagittal crest; increase in skull angularity; and enamel loss on the posterior side of molars 1 and 2 (Russell 1968b). Across all analyses, it is apparent that these subspecies are well-supported. Most subspecies remain distinct except for C. c. clarkii and C. c. parviceps. The type locality of C. c. clarkii is in Colorado (Table 1), so it is likely that samples included in this analysis are not representative of the genetic composition of the topotype, and these individuals have been influenced by admixture with other C. castanops groups in Texas (Fig. 4). Moreover, Bradley et al. (2023) advised eliminating the subspecific designation of ‘clarkii.’ Given that no other groups from Cratogeomys were investigated for this report, including sister groups C. goldmani or C. merriami, more cannot be said confidently regarding the taxonomy of the 7 subspecies found in Texas (Supplementary Data SD10). Further genomic investigation, either by genomic scans, other reduced representation techniques, or exome analysis is required that compare all groups within the genus (C. castanops, fumosus, goldmani, merriami, and planiceps; Hafner et al. 2004). However, Cytochrome b data (Bradley et al. 2023) show Cratogeomys species divided into 3 broad groups: (i) C. castanops and C. goldmani; (ii) C. fulvescens, C. merriami, and C. perotensis; and (iii) C. fumosus, C. gymnurus, and C. tylorhinus.

In a conservation management context, these data must be considered in light of what has already been published. Pocket gophers have a rich and substantial literature, but this is the first investigation into the population genomics of pocket gophers in Texas and adjacent New Mexico using large numbers of nuclear loci. This work suggests that a number of these groups are genetically similar despite phenotypic differences. And, given that only ~3% of any mammalian genome codes for protein (Mouse Genome Sequencing Consortium 2002), it is unlikely that many of the loci evaluated are functional and impact phenotype (Thibert-Plante and Hendry 2010). The genetic data provided are based on randomly selected, noncoding loci that should provide an unbiased estimate of genetic diversity (Hohenlohe et al. 2010). In other words, the methods used reveal overall genome-wide diversity but are likely to have missed coding variants that may be responsible for any morphological diversity that has been documented. Many groups display phenotypic plasticity, as between genera these groups are easily distinguished but within species they can be difficult to differentiate (Hall and Kelson 1959). For example, pleiotropy, when 1 gene influences other, often unrelated, phenotypic traits (Williams 1957) can influence phenotypic expression. A phenotypic trait thought to promote reproductive isolation would generate substantial differences in karyotype, and fixation of disruptive chromosomal rearrangements could lead to speciation (Davis et al. 1971; Patton and Yang 1977; Honeycutt and Schmidly 1979; Reig et al. 1989). This concern extends beyond pocket gophers to other subterranean groups such as Ctenomys (Ortells 1995), where chromosomal speciation has been speculated in a group that displays both morphological and allozymic similarity within its species.

In this first study using tens of thousands of genomic loci to investigate genetic similarity in nearly 40 groups in Texas, the following was determined based on the interpretation of the results and considering previous research: we show that reduced representation genome sequencing provides evidence that upholds current taxonomic classifications of Geomys species and C. castanops in Texas. In T. bottae, however, data show that classifications based on cranial morphology inflate the number of subspecies recognized, which could cause concerns in management decisions. Results support only 5 groups of T. bottae subspecies in Texas and neighboring New Mexico, typically viewed as having more than a dozen subspecies. Perhaps whole genome, transcriptome, or ectoparasite analyses could be conducted to strengthen the evidence of taxonomic relationships already established and offer possible avenues to clarify relationships that are not as resolved.

Such reductions were not necessary for the other 2 genera examined. Most Geomys species remain valid, apart from G. juggosicularis, which we combine with G. b. major. As mentioned in Sudman et al. (2006), G. breviceps in eastern Texas show interesting genetic results, and this group could benefit from targeted study with large sample sizes to determine if any of its currently recognized subspecies would qualify for elevation to species status. Finally, the only recommendation needed for C. castanops is to drop C. c. clarkii from the gopher taxonomic lexicon.

Despite the wealth of data from decades of research on these diverse subterranean rodents, our knowledge is not complete. Each taxonomic iteration provides researchers with confidence in groups that have remained relatively unchanged (C. castanops) or a “call to action” for those groups facing relentless restructuring (T. bottae). To conserve these groups properly, care must be taken regarding their classification in order to dissuade inaction of managing biodiversity.

Supplementary data

Supplementary data are available at Journal of Mammalogy online.

Supplementary Data SD1. BUSCO results for genome assemblies of G. bursarius and T. bottae.

Supplementary Data SD2. ddRADSeq sample sizes, locus counts, and private alleles for T. bottae.

Supplementary Data SD3. ddRADSeq sample sizes, locus counts, and private alleles for G. bursarius.

Supplementary Data SD4. ddRADSeq sample sizes, locus counts, and private alleles for Cratogeomys spp.

Supplementary Data SD5. Pairwise F_ST values between Thomomys bottae subspecies tested in analysis. Values greater than 0.25 are indicated in bold. The generic name and specific epithet have been dropped for readability. With only one representative, T. bottae pectoralis was not included.

Supplementary Data SD6. Pairwise F_ST values between Geomys groups tested in analysis. Because some of these are between true species, comparisons greater than 0.5 are indicated in bold.

Supplementary Data SD7. Pairwise F_ST values between Cratogeomys castanops subspecies tested in analysis. Values greater than 0.25 are indicated in bold. The generic name and specific epithet have been dropped for readability. With only a single representative, C. castanops perplanus was not included.

Supplementary Data SD8. Map of sampling localities for Thomomys spp.

Supplementary Data SD9. Map of sampling localities for Geomys spp.

Supplementary Data SD10. Map of sampling localities for Cratogeomys spp.

Acknowledgments

The completion of this project would not have been feasible without the help of numerous individuals, from graduate students and landowners to Texas Parks and Wildlife Department employees. Their contributions, especially those of all who contributed to work in the field, were substantial. We also thank the museum staff at the NSRL for quickly responding to our requests for museum specimens, particularly Heath Garner, Heidi Stevens, and Kathy McDonald. We also thank the sequencing team at Admera Health and BYU who produced the high-quality sequences that were ultimately the heart of this project. Also, we acknowledge the High-Performance Computing Cluster at Texas Tech University for use of their resources to facilitate analysis of these data. The authors thank Stony Brook Research Computing and Cyberinfrastructure, and the Institute for Advanced Computational Science at Stony Brook University for access to the high-performance SeaWulf computing system, which was made possible by a $1.4M National Science Foundation grant (#1531492).

Author contributions

MKH: conceptualization, data curation, formal analysis, investigation, methodology, validation, visualization, writing—original draft, writing—review & editing. EMR: data curation, formal investigation, methodology, resources. EAW: data curation, formal investigation, methodology, resources. TJS: data curation, formal investigation, methodology. LLL: data curation, formal investigation, methodology, resources. DMS: software, validation, visualization, writing—original draft, writing—review & editing. RMP: data curation, formal investigation, methodology, resources. LMD: software, validation, resources. RDB: conceptualization, data curation, funding acquisition, investigation, methodology, project administration, resources, supervision, validation, visualization, writing—review & editing. RDS: conceptualization, funding acquisition, project administration, resources, supervision, validation, writing—review & editing. DAR: conceptualization, data curation, funding acquisition, investigation, methodology, project administration, resources, supervision, validation, visualization, writing—review & editing.

Funding

State Wildlife Grant from the Texas Parks and Wildlife Department awarded to DAR, RDS, and RDB.

Conflict of interest

None declared.

Data availability

Data are available on the Sequence Read Archive (SRA) through the National Center for Biotechnology Information (NCBI) under the Project Number PRJNA1181766.

References

Appendix I

Specimens whose DNA was extracted for this study are listed below according to the NSRL Tissue and Karyotype number (TK) or collector code.

Specimens examined

Cratogeomys castanops angusticeps.—Texas; Maverick County, SSE on Loop 480 near Eagle Pass 28.66041831 N, 100.4842435 W, (TK197757, TK197752).

Cratogeomy castanops clarkii.—Texas; Brewster County, Elephant Mountain Wildlife Management Area, coordinates (TK197800). Texas, Presidio County, 35 miles S. Marfa 29.89241274 N, 104.0154314 W (TK197805), Texas; Hudspeth County, 1.5 miles S Sierra Blanca, 31.14690926 N, 105.3531065 W (TK197810). Texas; Pecos County 26.5 miles N Marathon, 30.472405, 102.9281231 (TK198763). Texas; Brewster County, Big Bend National Park, River Road East 29.18719751 N, 103.0146664 W (TK198765).

Cratogeomys castanops dalquesti.—Texas; Sterling County, Highway 163 W of Sterling City, 31.953721 N, 101.035489 W (TK198772). Texas; Sterling County, Highway 163 W of Sterling City, 31.876122, 101.043423 W (TK198773). Texas; Sterling County, Highway 163 W of Sterling City near Willow Creek Bridge, 31.966355 N, 101.022844 W (TK198774).

Cratogeomys castanops lacrimalis.—New Mexico; Chaves County, Roswell along railroad Track, 33.37764 N, 104.50393 W (TK197951). New Mexico; Chaves County, Roswell along railroad Track, 33.36185 N, 104.49682 W (TK197952). New Mexico; Chaves County, Roswell along railroad Track 33.36922 N, 104.49770 W (TK197953).

Cratogeomys castanops parviceps.—New Mexico; Otero County, Fort Bliss, 32.53488 N, 105.93144 W (TK51761). New Mexico; Otero County, Fort Bliss, no coordinates reported (TK51805). New Mexico; Otero County, 18 miles SW Alamogordo, 32.74229 N, 106.19960 W (TK197968). New Mexico; Otero County, 12 miles SW Alamogordo, 32.77805 N, 106.17263 W (TK197973).

Cratogeomys castanops perplanus.—Texas; Oldham County, Boys Ranch, 35.52892 N, 102.26096 W (TK198778).

Cratogeomys castanops tamaulipensis.—MEXICO: Tamaulipas; 0.86 miles upriver from old bridge on Rio Grande in Matamoros, no coordinates reported (TK27005, TK27006, TK27007).

Geomys arenarius arenarius.—Texas; El Paso County. 1 mile S, 0.25 miles W Fabens, 31.65853 N, 108.37555 W (TK47533, TK47534, TK47536, TK47537, TK47538, TK47540, TK47541).

Geomys arenarius brevirostris.—New Mexico; White Sands National Monument, 32.79117 N, 106.22664 W (from A.N. Kozora 2011 M.S. thesis ANH2, ANH3, ANH4, ANH9).

Geomys attwateri.—Texas; Aransas County, 1.1 miles E Rockport 28.0313 N, 97.05155 W (TK200418, TK200419, TK200420). Texas; Aransas County, 0.57 miles SE Rockport 28.02743333 N, 97.04523333 W (TK200421, TK200422, TK400423).

Geomys bursarius ammophilus.—Texas; Dewitt County, 1.13 miles S Cuero 29.08288333 N, 97.27626667 W (TK200269, TK200270, TK200271, TK200272).

Geomys bursarius brazensis.—Texas; Grimes County, 5.81 miles E Kurten on CR 101, 30.8151 N, 6.1708 W (TK200288, TK200289, TK200290, TK200291, TK200292, TK200293, TK200443). Texas; Grimes County, 7.99 miles, SE Kurten on CR 162, 30.7159 N, 96.1549 W (TK200444).

Geomys bursarius major.—Texas; Donley County, TX-70 S of Clarendon, 34.876098 N, 100.877093 W (TK1978755). Texas; Donley County, TX-70 S of Clarendon, 34.876122 N, 100.876892 W (TK1978756), Texas; Donley County, W of Clarendon on FM 2362, 34.947064 N, 100.974942 W (TK1978757).

Geomys bursarius dutcheri.—Oklahoma, Muskogee County, Fort Gibson Historic Site, 35.8043 N, 95.2578 W (TK200185, TK200186, TK200187, TK200188, TK200189).

Geomys breviceps sagittalis.—Texas; Milam County, 1.5 miles W Maysfield, 30.92020 N, 96.84376 W (TK24932), Texas; Wood County, 3.5 miles SE Quitman (TK48934). Texas; Wood County 2 miles SE Quitman 28.90088, 99.05751 (TK48935, TK48940). Texas; Wood County 2 miles S Quitman, 32.79324 N, 95.40955 W (TK48939).

Geomys juggosicularis.—Texas; Dallam County, 12 miles NE Texline, Thompson Grove Campground, Rita Blanca National Grassland 36.41508 N, 102.80576 W (TK197190, TK197202). Texas; Dallam County, 12 miles NE Texline, Thompson Grove Campground, Rita Blanca National Grassland, 36.41501 N, 102.86021 W (TK197204). Texas; Dallam County, 12 miles NE Texline, Thompson Grove Campground, Rita Blanca National Grassland, 36.41454 N, 102.84687 W (TK197201, TK197205, TK197206).

Geomys knoxjonesi.—Texas; Yoakum County, Yoakum Dunes Wildlife Management Area 33.40539 N, 102.69477 W (TK197893). Texas; Yoakum County, Yoakum Dunes Wildlife Management Area, 33.40362 N, 102.70198 W (TK197894). Texas; Winkler County, 4.5 miles N, 4.5 miles E Kermit, 31.93549 N, 103.03047 (TK200769, TK200771). Texas; Winkler County, 5.1 miles N, 5.1 miles E Kermit, 31.94448 N, 103.01316 W (TK200770)

Geomys personatus davisi.—Texas; Zapata County. Zapata Airfield 26.964846 N, 99.249259 W (TK197860, TK197861, TK197862, TK197863).

Geomys personatus maritimus.—Texas; Nueces County, Corner of Flour Bluff Road and Graham Rd 27.82884 N, 97.10864 W (TK28648, TK28649). Texas; Nueces County, Flour Bluff, Corpus Christi 27.64219 N, 97.31947 (TK197942).

Geomys personatus megapotamus.—Texas; 4.5 miles SE Oilton 27.457124 N, 98.922195 W (TK197753, TK197755, TK197756, TK197757, TK197758, TK197759, TK197760 TK197761).

Geomys personatus personatus.—Texas; Nueces County, Mustang Island, 27.73149 N, 97.09115 W (TK943920, TK943922, TK943927).

Geomys streckeri.—Texas; Dimmitt County, Carrizo Springs, 33.628626 N, 101.756913 W (TK197851, TK197852, TK197853, TK197854, TK197855, TK197858, TK197858).

Geomys texensis bakeri.—Texas; Uvalde County, 11 miles S Sabinal near Sabinal River on FM 187 29.153836 N, 99.480423 W (TK197763, TK197764, TK197765, TK197766, TK197768, TK197769, TK197770).

Geomys texensis llanensis.—Texas; Llano County, Llano River Golf Course 2.11 miles W Llano on FM 152, 30.74798333 N, 98.70911667 W (TK200139). Texas; Medina County, 0.97 miles E D’Hanis 29.3318166 N, 99.26471667 W (TK200341, TK200342, TK200343, TK200344, TK200345).

Geomys texensis texensis.—Texas; Mason County, Mason Mountain Wildlife Management Area, 30.82652 N, 99.21765 W (TK160809, TK160810, TK111633, TK111636, TK111639, TK133623, TK136592, TK116745, TK147179).

Geomys tropicalis.—MEXICO: Tamaulipas; 2.17 miles S Altamira, no coordinates reported (TK27097, TK27107, TK27108, TK27115, TK27116, TK27117, TK27118, TK27119).

Thomomys bottae actuosus.—New Mexico; Lincoln County, Corona, 34.23691 N, 105.59778 W (TK200930, TK200935).

Thomomys bottae confinalis.—Texas; Kimble County, 1.13 miles SW London, 30.62819 N, 99.59002 W (TK185684). Texas; Kimble County, 1.75 miles S London, 30.65603 N, 99.58696 W (TK185685), Texas; Kimble County 2.25 miles SSW London, 30.64806 N, 99.58362 W (TK26872, TK26873). Texas; Kimble County 2.5 miles SSW London, 30.64806 N, 99.58362 W, (TK26874). Texas; Kimble County, 3 miles SSW London, 30.64127 N, 99.58683 W (TK26995). Texas; Kimble County, 4 miles SSW London 30.62768 N, 99.59321 W (TK26996). Texas; Kimble County 2 miles S London, 30.65255 N, 99.56746 W (TK27171). Texas; Kimble County, 1.5 miles S London, 30.65992 N, 99.56750 W (TK27172). Texas; Crockett County, 14 miles N, 16 miles W, Ozona, 30.90874 N, 101.48440 (TK26903).

Thomomys bottae guadalupensis.—Texas; Culberson County. Guadalupe Mountains National Park, McKittrick Canyon, 31.978057N, 104.753360 W (TK200912, TK200913, TK200914). Texas; Culberson County. Guadalupe Mountains National Park, McKittrick Canyon Trail, 31.893659 N, 104.825389 W (TK197811).

Thomomys bottae lachuguilla.—Texas; Franklin Mountains State Park, Tom Mays Unit, 31.928609 N, 106.508308 W (TK198760, TK198761)

Thomomys bottae limpiae.—Texas; Jeff Davis County, 10.7 miles NE Fort Davis, no coordinates reported (TK75201). Texas; Jeff Davis County, 10.9 miles NE Fort Davis, no coordinates reported (TK75209). Texas; Jeff Davis County, 11.3 miles NE Fort Davis, no coordinates reported (TK75205, TK75207). Texas; Jeff Davis County, 13.8 miles NE Fort Davis, no coordinates reported (TK75201, TK75206). Texas; Jeff Davis County, 15.8 miles NE Fort Davis, no coordinates reported (TK75208). Texas; Jeff Davis County, 1 mile S Wild Rose Pass on Texas State Highway 17, 30.69012 N, 103.78947 W (TK197808).

Thomomys bottae limitarus.—Prior to Beauchamp-Martin et al. (2019), the following specimens were classified as T. b. limitarus. All specimens used in the analysis have been recategorized and its new classification is listed in brackets after the TK number.

Texas; Jeff Davis County, Mt. Livermore Preserve Park, no coordinates provided (TK84861 [T. b. texensis]). Texas; Brewster County, Big Bend Ranch State Park 29.44233 N, 103.78142 W (TK54878, TK54879, TK54880 [T. b. pervarius]). Texas; Upton County, 1.5 miles E McCamey, 31.13335 N, 102.19451 W (TK26990 [T. b. robertbakeri]). Texas; Upton County, 4 miles N, 4 miles E McCamey, 31.19326 N, 102.15371 W, (TK26992 [T. b. robertbakeri]). Texas; Upton County, 3 mi E McCamey, 31.72650 N, 103.51314 W (TK26807 [T. b. robertbakeri]). Texas; Reagan County, 1 mi W Best, 31.23397 N, 101.61865 (TK26902 [T. b. robertbakeri]). Texas; Reagan County, 3 mi W Big Lake, 31.19210 N, 101.50974 W (TK27169 [T. b. robertbakeri]). Texas; Irion County, 0.5 mi W Barnhart, 31.14461 N, 101.16218 W (TK27170 [T. b. robertbakeri]). Texas, Crockett County, 14 miles N, 16 miles W Ozona, 30.90874 N, 101.48440 N (TK26903 [T. b. robertbakeri]). Texas; Jeff Davis County, Mt. Livermore Preserve Park, no coordinates provided (TK84861 [T. b. texensis]). Texas; Brewster County, Big Bend Ranch State Park 29.44233 N, 103.78142 W (TK54878 [T. b. pervarius]).

Thomomys bottae pectoralis.—New Mexico; Eddy County, 7 miles E of Queen at Junction of NM 147 and CR 410, 32.20234 N, 104.62622 W (TK200917).

Thomomys bottae pervarius.—Texas; Presidio County. Big Bend State Natural Area, no coordinates reported.

Thomomys bottae robertbakeri.—Texas; Upton County 2.5 miles E McCamey, 31.13377 N, 102.17741 W (TK26877). Texas; Upton County, 1.5 miles E McCamey, 31.13335 N, 102.19451 W (TK26990, TK26991). Texas; Upton County, 12 miles N, 5 miles E McCamey, 31.31137 N, 102.14048 W (TK27166), Texas; Upton County, McCamey County Club, 31.13271 N, 102.22016 W (TK27168).

Thomomys bottae ruidosae.—New Mexico; Otero County, Lincoln County National Forest, Scott Able Rd, 32.746699 N, 105.708907 W (TK198775). New Mexico; Otero County, Lincoln County National Forest, 32.763391 N, 105.72574 W (TK198776). New Mexico; Otero County, Lincoln County, Ruidoso, White Mountain Elementary School Athletic Fields, 33.35695 N, 105.66016 W (TK200926). New Mexico; Otero County, Lincoln County, Ruidoso, Schoolhouse Park, 33.32211 N, 105.64202 W (TK200929, TK200932, TK200933, TK200934). New Mexico; Otero County, Fort Bliss, 32.47839 N, 105.74768 W (TK51802). New Mexico; Otero County, Lincoln County, Lincoln County National Forest, Deerhead Campground (TK49857, TK49858, TK49859, TK49860), no coordinates reported.

Thomomys bottae scotophilus.—Texas; Brewster County, Sierra Diablo Wildlife Management Area, 31.437333 N, 104.908389 W (TK197878, TK198777). Texas; Brewster County, Sierra Diablo Wildlife Management Area, 31.26674 N, 104.89892 W (TK200891, TK200894, TK200895). Texas; Brewster County, Sierra Diablo Wildlife Management Area, 31.27056 N, 104.90970 W (TK54174, TK54190, TK54210).

Thomomys bottae spatiosus.—Texas; Brewster County, Elephant Mountain Wildlife Management Area, 30.037050 N, 103.540591 W (TK199788, TK199810, TK199811, TK199831, TK200855). Texas; Pecos County 22 miles N of Marathon, 30.40809 N, 102.984812 (TK198764).

Thomomys bottae tularosae.—New Mexico; 6.8 miles S Cloudcroft, 0.7 miles W on Park Rd 634, 32.89685 N, 105.78082 W (TK51810). New Mexico; 6.8 miles S Cloudcroft, 0.5 miles W on Park Rd 634, 32.86466 N, 105.79669 W (TK51811).

Thomomys bottae texensis.—Texas; Jeff Davis, Mt. Livermore Preserve. 30.66372 N, 104.15809 W (TK84858, TK84859, TK84861, TK83553, TK83554, TK83555, TK83589, TK78999, TK79072, TK79072, TK79075). Texas; Jeff Davis, Davis Mountains State Park, 0.85 miles W show tank, 30.67113 N, 104.14266 W (TK92798). Texas; Jeff Davis, Mt. Livermore Preserve, road towards cattle tank, 30.65736 N, 104.16654 W (TK79003).