Microbial communities are diverse assemblages of microbiotic species that, through interactions with each other and with the physical and chemical components of their abiotic environments, have substantial impacts on global processes. Microbes play an important role in global nutrient cycling (Treseder et al., 2016; Heijboer et al., 2018) and energy flow through ecosystems (Konopka, 2009). In turn, microbial communities are structured by the physical and chemical properties (Leff et al., 2015; Garcia et al., 2020; Xia et al., 2020) of the soil substrate, including soil moisture, C:N ratio, pH, and total carbon content (Shen et al., 2013; Li et al., 2017).

In addition to their interactions with abiotic processes, soil microbiota structure biotic diversity and regulate the health of hosts that house the microbial communities (Ichinohe et al., 2011; Shen et al., 2018). Soil microbes influence plant and animal communities (Lau and Lennon, 2011; Seastedt et al., submitted¹) through mechanisms such as increasing plant nutrient acquisition (Hestrin et al., 2019) and resistance to desiccation (Xi et al., 2018) and inhibiting or facilitating the establishment of pathogens (Perez et al., 2008; van Elsas et al., 2012). Soil microbes are in turn governed by the actions of plants (Zak et al., 2003; Prescott and Grayston, 2013; Lange et al., 2015) and animals (Kandeler et al., 1999; Cline et al., 2017; Bray et al., 2019), creating feedbacks between soil microbial and aboveground communities (Bartelt-Ryser et al., 2005).

Biotic and abiotic processes that influence soil characteristics may be predicted to govern microbial diversity. For instance, ecosystem engineers that influence sediment abiotic properties (e.g., bioturbating shrimp, Laverock et al., 2010; Populus, Ciadamidaro et al., 2013) or soil nutrients (e.g., prairie dogs, Anacker et al., 2021) should thus also determine the microbial communities present (Gutiérrez and Jones, 2006; Cregger et al., 2018; Zotti et al., 2020). Similarly, mass mortality events in animals supply nutrient pulses that should alter microbial communities and contribute to terrestrial nutrient cycling (Metcalf et al., 2016b). Mass mortality in ecosystem engineers or keystone species, which influence the abundance of other (typically plant and animal) taxa, could have an especially pronounced effect. Soil microbiota can regulate the microbial pathogens causing such mass mortality, for instance if soil microbial communities contain animal pathogens or reservoirs for animal pathogens (Markman et al., 2018) or, conversely, microbes that inhibit establishment of animal pathogens. Through facilitation or inhibition of pathogens (Perez et al., 2008; van Elsas et al., 2012), soil microbes thus contribute to the maintenance of biodiversity of plants and animals.

Black-tailed prairie dogs (Cynomys ludovicianus) are social, fossorial ground squirrels inhabiting North American grasslands that build extensive underground burrows. Burrows typically range in length from 5 to 10 m long and extend as deep as 3–4 m below ground (Wilcomb, 1954; Hoogland, 1995). Burrows maximize air and water flow through the burrow and minimize water retention within the burrow, thus creating a moist but not wet environment. Their burrows increase soil porosity (Gedeon et al., 2012, which can facilitate deeper penetration of precipitation (Munn, 1993). Prairie dogs also increase the total nitrogen content and productivity of soils inside or near their burrows, leading to higher plant growth and diversity (Whicker and Detling, 1988; Holland and Detling, 1990).

More than half of a prairie dog’s life is spent within its burrow: prairie dogs use their burrows for reproducing, storing food, and escaping from both predators and the environment (Hoogland, 1995). Therefore, burrows are complex and heterogeneous in structure, and include spatially segregated chambers with various purposes, including nesting, hibernating (in species that hibernate; Cooke and Swiecki, 1992), defecating, and burying or isolating dead kin (Burns et al., 1989). Prairie dogs can die in their burrows over winter as a result of insufficient resources, and at other times of year from causes such as infectious disease. The primary disease affecting prairie dogs is sylvatic plague, caused by the Gram-negative bacterium Yersinia pestis. Typically transmitted by fleas, the pathogen is extremely virulent to prairie dogs, with individual colonies undergoing severe population declines ranging between 85% and complete extinction (Cully et al., 2010). These die-offs can thus result in hundreds of kilograms of carcasses appearing over the course of several weeks. In between epizootics, the plague reservoir is unknown: Some have hypothesized the pathogen persists in an alternative mammalian (Salkeld et al., 2010) host or flea vector (Webb et al., 2006) while others have posited that the reservoir is telluric (Drancourt and earlier; Eisen et al., 2008), residing in an invertebrate such as a nematode or amoeba (Markman et al., 2018).

Prairie dogs regularly clean out their burrows, leaving piles of nesting material, scat, and bones near some entrances of burrows. This excavated soil provides an opportunity to non-invasively explore the microbial composition of various locations within prairie dog burrows. We hypothesize that prairie dogs structure soil microbial communities through their functional partitioning of burrows, and that this structure may be pronounced after mass mortality caused by the pathogen Y. pestis. This study is the first to characterize the fine-scale spatial variation in microbial communities in the complex structure of prairie dog burrows.

We obtained 79 soil samples from 64 burrows in 6 colonies. Nitrogen content averaged 0.267% (range 0.062–0.685%) and carbon content averaged 3.56% (range 0.645–7.59%); mean C:N ratio was 14.8% (range 9.18–41.5%). Mean soil water content was 0.075 g/g (range 0.004–0.22) and mean pH was 7.89 (range 6.18–9.06). There was significant variation in soil properties among colonies (Supplementary Figure 1) and among burrow regions (Supplementary Figure 2). In particular, soil moisture was significantly higher in colony 30A and lower in colony 1A than other sites, and pH was significantly lower in colony 12A than several other sites (but sample sizes in 12A and 30A were small). Soil moisture was significantly higher in soil collected from the burrow entrance than in excavated soil containing prairie dog bones. The C:N ratio was significantly higher in soil sampled from excavated soil containing bones and scat than all other regions except those with scat. Total carbon, total nitrogen, and pH did not vary across sampling regions.

All clustering methods produced highly similar results, with the UniFrac unweighted method resulting in the highest proportion of variance explained by the first three axes. Samples collected from recently excavated soil adjacent to burrows clustered slightly on Axis 1, but samples from different regions were largely overlapping (Supplementary Figure 3).

The best initial model (excluding single taxa) of Shannon’s ecological diversity included the predictors: burrow region, pH, water content, and interactions between pH and water content and between carbon and nitrogen content (AIC 254.98). Colonies did not differ in ecological diversity, and inclusion of colony as a predictor worsened the model (AIC 263.65). Inclusion of Yersinia presence as a predictor worsened the model, but not significantly (AIC 256.97). All variables in the model significantly influenced diversity (Supplementary Table 1). Diversity was lowest in soil collected in the presence of scat, followed by soil with bones and scat, and was highest in soil recently excavated from burrows and from those with plague-positive animals (Figure 2). The best model of phylogenetic diversity included the same predictor variables (in this case, inclusion of colony as a predictor led to a worse, but statistically indistinguishable model: AIC with colony = 1125.9, AIC without colony = 1125.7). Inclusion of Yersinia presence as a predictor resulted in a statistically indistinguishable model (AIC 1125.8). All predictors significantly influenced phylogenetic diversity except for burrow location, which had a marginally significant effect (p = 0.074). Similar to the pattern observed for ecological diversity, phylogenetic diversity exhibited a trend toward lower diversity in soil collected in the presence of scat, followed by soil with bones and scat, and higher diversity in soil recently excavated from burrows, soil from burrows inferred to contain dead animals, and soil excavated from burrows with plague-positive animals (Figure 2).

Adding the relative abundance of Enterobacteriaceae improved the AIC of both models (Shannon diversity AIC 252.05, significant improvement; Faith diversity AIC 1124.3, marginal improvement). The relative abundance of Enterobacteriaceae had a negative effect on both ecological and phylogenetic diversity, although this effect seemed to be driven by an outlier with a relatively high proportion of Enterobacteriaceae and low diversity. Removing the outlier changed the magnitude (and significance) of the relationship, but the trend toward an inverse relationship persisted. The improvement of model fit with the inclusion of relative abundance of Enterobacteriaceae was not unique to this family; in fact, the inclusion of 157 single taxa significantly improved model fit (reducing AIC by more than 2); 59 families (one at a time) reduced AIC values by >10. In particular, the best single-taxon model of Shannon’s diversity included relative abundance of Planococcaceae (AIC 188.53) in addition to the previous predictors, and the only other model within 10 AIC was a model including relative abundance of Micrococcaceae (AIC 194.10). Both of these taxa exhibited a strong negative relationship with Shannon’s diversity. Similarly, the inclusion of 157 single taxa significantly improved model fit (reducing AIC by more than 2) for phylogenetic diversity, and 67 families reduced AIC values by >10. The best single-taxon model of phylogenetic diversity included the relative abundance of an unknown family in order WD2101 (class Phycisphaerae, AIC 1064.83) in addition to the previous predictors, and the only other model within 10 AIC was a model including relative abundance of an unknown family in order iii1–15 (class Acidobacteria-6, AIC 1071.36). Both of these taxa exhibited a positive relationship with phylogenetic diversity.

Burrow regions differed significantly in the most abundant taxa at all taxonomic levels (phylum, class, order, family, and genus; Supplementary Figure 4). Across all samples, the dominant family averaged 12.5% of the total sequences per sample, and ranged from comprising 5.6–49.2% of the total sequences per sample. In soils collected from burrows inferred to currently contain dead animals, Firmicutes were more abundant than expected (Chi-square = 74.528, df = 30, p-value = 1.172e−05). Burrows with dead animals contained more Bacilli (and Bacillales) and Rubrobacteria (and Rubrobacterales) than expected, while soils containing bones were characterized by a lower abundance of Alphaproteobacteria than expected (Chi-square = 182.36, df = 72, p-value = 1.545e−11). Soils containing bones and scat possessed a lower abundance of Rhizobiales than expected (Chi-square = 234.44, df = 96, p-value = 1.382e−13).

Forty-eight bacterial families and 76 bacterial genera varied significantly in relative abundance across burrow regions (Supplementary Tables 2–5). Among the taxa most significantly varying across burrow regions were an unknown family and genus in the AKIW781 order of class Chloroflexi, which was an order of magnitude higher in soil with bones and scat (Figure 3A); Deinococcus (Deinococcaceae), which was an order of magnitude higher in soil with bones and scat and an order of magnitude lower in soil associated with dead animals (Figure 3B); an unknown genus in Planococcaceae, which was highest in soil associated with dead animals (Figure 3C); and Cellulosimicrobium (Promicromonosporaceae, Actinomycetales), which was highest in soils sampled with scat (Figure 3D).

Enterobacteriaceae, the family containing Y. pestis, was found in all soil samples, but at low proportions (never exceeding 3%). The proportion of Enterobacteriaceae sequences was significantly higher in samples with higher C:N (p = 0.0002) and in burrow regions associated with dead animals than in other regions (p = 0.013). Although we ran this model first due to our particular interest in the family, we also aimed to determine the extent to which spatial variation in abundance was characteristic shared by many microbial taxa. When we ran separate models for all 565 families and 990 genera, the false discovery rate correction led to a loss of statistical significance for spatial variation in Enterobacteriaceae (data not shown). Yersinia was identified in four samples from two burrows in one colony (19A). All Yersinia-containing samples were collected from waste chambers. Presence of this genus was a significant predictor of the relative abundance of five bacterial families and eight microbial genera (Figure 4 and Tables 1, 2). All of these taxa were found in significantly higher abundance in samples where Yersinia was present. Many of these genera (e.g., 9 out of the 10 strongest associations) were extremely rare taxa that appeared only or primarily in the samples containing Yersinia. The presence of Yersinia in a sample was associated with slightly, but not significantly, lower microbial diversity within samples (Shannon without Yersinia 8.303, Shannon with Yersinia 8.059, p = 0.14; Faith PD without Yersinia 81.807, Faith PD with Yersinia 75.111, p = 0.17; Figure 5).

Microbial communities as a whole varied – and many specific taxa differed in relative abundance – at small spatial scales among regions of a prairie dog burrow following a mass mortality event. Dominant taxa were consistent with predictions of microbial succession following the nutrient pulse that occurs during decomposition of mammalian corpses (Metcalf et al., 2016a, b). In addition, several taxa were significantly associated with the presence of Yersinia in soil samples, primarily as a result of taxa of low abundance found at higher abundance when Yersinia was present. Both ecological and phylogenetic diversity resulted from the combined influences of soil properties and burrow region.

Other studies have shown similar degrees of fine-scale spatial structure in microbial communities resulting from niche differentiation (Zhuang et al., 2020), particularly in microbial communities associated with plant roots (Aas et al., 2019) and other plant tissues (Cregger et al., 2018). Niche diversification may be particularly likely when niches are divergent even at small spatial scales, when specific microbes present in high abundance in certain environments exert selection on other microbial taxa (e.g., predatory microbes) or when microenvironments are less hospitable (e.g., very dry). In this system, microbial communities associated with scat may be specialized for living in the mammalian gut, metabolizing plant tissues, or both. Soil collected with bones were the driest soils we sampled, thereby potentially exerting strong selection on microbial communities in these soils.

Fine-scale spatial structure could also arise from community assembly (Nemergut et al., 2013) and succession processes such as colonization of a deceased animal from soil microbiota (Metcalf et al., 2016b), particularly if animals died in a spatially structured way or were moved to specific locations after death – scenarios that are consistent with the few existing observations of deceased prairie dogs within burrows (Burns et al., 1989). For instance, Enterobacteriaceae are abundant in the early stages of corpse decay, while Planococcaceae become more abundant as corpse decay progresses (Metcalf et al., 2016a). This is consistent with our observation of significantly higher abundance of both taxa in soils collected near dead animals.

Keystone microbial taxa (Banerjee et al., 2018) can influence the abundance of other community members based on ecological interactions (Herren and McMahon, 2018) including the prevention of pathogen establishment (Trivedi et al., 2017). We found >50 taxa that significantly influenced ecological or phylogenetic diversity among samples, with some having particularly strong effects. Four single taxa [Planococcaceae, Micrococcaceae, unknown family in WD2101 (Planctomycetes), and unknown family in iii1–15 (Acidobacteria)] were statistically separated as predictors of diversity (in conjunction with abiotic soil properties) from other taxa, indicating their potential role as keystone taxa. A negative relationship between Planococcaceae and ecological diversity supports previous findings of this family becoming more abundant after disturbance of an ecological community (Aanderud et al., 2019). Micrococcaceae have been associated with increased plant growth (Hong et al., 2016), which could cause feedbacks with microbial diversity, although the mechanism underlying this potential relationship is not clear. Both WD2101 and iii1–15 are among the most abundant soil bacteria globally (Delgado-Baquerizo et al., 2018). The lack of taxon sharing within WD2101 (<2% shared OTUs) even in similar environments (Dedysh et al., 2020) and the low degree of genomic match to characterized sequences (Delgado-Baquerizo et al., 2018) suggest a large amount of cryptic diversity in the group that could be a driving force behind the high phylogenetic diversity we found here. The abundance of iii1–15 responds to soil moisture (Barnard et al., 2013), which could provide a mechanism for its relationship with phylogenetic diversity (Brockett et al., 2012).

Among the taxa that varied spatially within prairie dog burrows was an unknown member of the AKIW781 class (order Chloroflexi), found here with bones and scat, which has previously been described in soils from deserts in North and South America (Mogul et al., 2017; Lucas et al., 2020) and is likely adapted to dry conditions. Similarly, we found Deinococcus to be higher in soils with bones and scat, which may be not only drier but more exposed to sunlight than soils excavated from other parts of the burrow. Deinococcus is resistant to solar radiation and increases in relative abundance in irradiated soils (Ogwu et al., 2019). An unknown genus in Planococcaceae was highest in soils associated with dead animals, consistent with previous description of the abundance of this family in later stages of the decomposition process (Metcalf et al., 2016a). Finally, Cellulosimicrobium was found at highest relative abundance in soils containing scat, which supports the role of this genus in breaking down plant material (Bakalidou et al., 2002; Schumann and Stackebrandt, 2015).

In line with other studies of pathogens and soil microbial diversity (van Elsas et al., 2012), the presence of Yersinia in a sample was negatively associated with microbial diversity (although the relationship here was not significant). The most notable microbial association with Yersinia was with Dictyostelium, an amoeba that consumes bacteria. Previous experimental work has shown that Y. pestis can escape phagocytosis by and replicate within D. discoideum for at least 48 h (in comparison with control bacteria, which were consumed within 1 h; Markman et al., 2018). The prevalence of Dictyostelium (present in 2 of 158 samples) and another amoeba, Acanthamoeba (10 of 158 samples), in our soils was lower than that recovered in Markman et al. (2018), although the methods of recovery differed. To our knowledge, this is the first study to detect Yersinia in soil samples collected from prairie dog colonies. Although we were unable to classify the sequences at the species level due to read length constraints, this suggestive finding adds to the collective evidence that Y. pestis is present in prairie dog colonies in the absence of epizootics (3 years after the prairie-dog population die-off) and that soil amoebae may be a potential reservoir for plague in inter-epizootic intervals.

Our results show that variation in soil microbial communities occurs at fine spatial scales in relation to functional partitioning of below-ground space by a social mammalian herbivore. This fine-scale structure likely interacts with mass mortality events, for example by sudden drastic increases in input to certain physical burrow regions (e.g., chambers used for quarantining dead individuals). The existence of fine-scale spatial structure in community diversity in this and other studies suggests that estimates of beta-diversity should account for fine-scale structure in order to accurately estimate the true degree of diversity. Collectively, our results demonstrate how soil microbial communities can interact with animal pathogens (van Elsas et al., 2012; Trivedi et al., 2017) to shape above- and below-ground biodiversity in grasslands.

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://qiita.microbio.me/emp, 11519; https://www.ebi.ac.uk/ena/browser/view/ERP106314.

LC-S conceived, designed, and oversaw this study. CK and LC-S analyzed the data and wrote the manuscript. Both authors contributed to the article and approved the submitted version.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.