The raw sequences of 70 fecal samples totaled 5,407,660, with an average of 77,252 reads per sample (range 28,178 - 357,624; SD = 41,072). After quality filtering (Q > 30) and denoising, resulting on a grand total of 3,939,947 non-chimeric sequences, averaging 56,285 reads per sample (range 15,839 – 256,779; SD = 30,147). We detected a total number of 5,118 ZOTUs (zero-radius operational taxonomic units) by clustering all these reads based on a 100% similarity threshold. Observed OTUs leveled off at a sequencing depth of 40,000 indicated by the fact that enough OTUs have been detected to adequately characterize the microbial communities and the number of reads is not a limiting factor for OTU detection beyond 40,000 reads (Supplementary Information: Fig. S1). Some samples were characterized by a higher percentage of bacterial species than others and a shallow gradient from all samples indicated the relative abundance and incidence of those species were more evenly distributed. This observation was supported by the species rank abundance plot and the incidence abundance plot, as they showed a similar pattern across the dataset (Supplementary Information: Fig. S2A-2B).
Fecal microbiota composition
Overall, 11 phyla, 57 families, 114 genera, and 220 species were found in all samples of Cynomys ludovicianus. Bacteroidetes was the most abundant phyla representing 57.3% of the fecal samples followed by Firmicutes (39.6%), Proteobacteria (1.5%), Tenericutes (0.9%), and Actinobacteria (0.4%) (Fig. 2A). The top ten families, Rikenellaceae (19.2%), Bacteroidaceae (16.7%), Prevotellaceae (13.4%), Unclassified_Clostridiales (12.5%), Clostridiaceae (10.6%), Eubacteriaceae (7.3%), Porphyromonadaceae (6.5%), Lachnospiraceae (6.3%), Lactobacillaceae (1.3%), and Ruminococcaceae (1.2%), made up 95% of the fecal microbiota (Fig. 2B). The most abundant genera dominated in all samples including Alistipes (19.2%), Bacteroides (16.5%), Prevotella (12.2%), Eubacterium (7.3%), Clostridium (6.4%), Blautia (5.8%), Parabacteroides (4.2%), Anaerovorax (3.8%), Roseburia (3.3%), and Alkaliphilus (1.7%). Microbial communities were dominant by five species, Alistipes shahii (13.1%), Prevotella shahii (10.4%), Bacteroides rodentium (7.3%), Eubacterium oxidoreducens (3.9%), and Anaerovorax odorimutans (3.8%) with the remaining species accounting for less than 5% of the mean relative abundance.
Fecal microbiota alpha diversity
The number of observed OTUs varied from 647 to 1922 per sample (Supplementary Information: Table S1). A strong linear relationship was found between the observed OTU richness and the total read counts of the sample in the dataset (Welch Two Sample t-test: t = 14.349, df = 69.077, p < 0.001; Fig. 3). Microbial alpha diversity was summarized as observed OTU richness (R = 1474.3±25.54), Shannon diversity index (H’ = 6.45±0.04), as well as Faith’s phylogenetic diversity (PD = 11.55±0.14) which was followed by the analysis of variance among sampling sites. There were no significant spatial autocorrelations in Moran's I test for alpha diversity metrics (Observed: p = 0.838; Shannon: p = 781; Faith’s PD: p = 0.522). The distribution of observed OTUs was significantly varied across sites (ANOVA, F = 9.90, p < 0.01; Fig. 4A, Table 2). Significant differences were present between alpha diversity estimates of five different experimental groups (Shannon: F = 5.27, p = 0.00097; Faith’s PD: F = 5.01, p = 0.0013; Fig. 4B and 4C). There was no statistical difference between Dallam and Randall regarding intra-individual diversity (Fig. 4A-4C). For all three alpha diversity metrics, diversity and richness were significantly greater in Dallam, while lower in Bailey (Tukey’s HSD test, p<0.05; Fig. 4A-4C). We found that adding habitat as a covariate in the linear regression model increased the association between alpha diversity and sampling sites (8.0%), indicating part of the association was explained by habitat (p < 0.001; Table 3). There was a evidence of a multiple regression relationship between alpha diversity and the two predictor variables (F(4,65) = 9.901, R2 = 0.38, p < 0.001). Significant variation was observed in alpha diversity between habitats (Observed: F = 31.38, p < 0.001; Shannon: F = 14.50, p = 0.0003; Faith’s PD: F = 12.45, p = 0.0007; Table 3). Diversity was significantly higher in rural areas compared to their counterparts in urban areas (Tukey’s HSD test, p < 0.001; Fig. 5A-5C).
Table 1 Sample information including bioclimatic data
Sampling sites
|
Habitat
|
Elevation (m)
|
Cumulative precipitation (cm)
|
Average minimum temperature (°F)
|
Average maximum temperature (°F)
|
Lubbock (n=12)
|
Urban
|
970
|
19.1
|
63.6
|
91.5
|
Hockley (n=14)
|
Urban
|
1010
|
19.8
|
62.8
|
91.2
|
Bailey (n=16)
|
Urban
|
1163
|
18.4
|
61.4
|
90.3
|
Randall (n=15)
|
Rural
|
1107
|
20.1
|
62.7
|
89.9
|
Dallam (n=13)
|
Rural
|
1320
|
22.2
|
59.9
|
88.5
|
Table 2 Alpha diversity estimates in fecal microbiome of Cynomys ludovicianus across sampling sites and habitat
Sampling sites
|
Samples
|
Total reads
|
Reads after filtering
|
Observed OTU’s
|
Shannon
|
Faith’s PD
|
Lubbock
|
n=12
|
69949±5780
|
51434±4548
|
1481.1±116.8
|
6.49±0.15
|
11.79±0.462
|
Hockley
|
n=14
|
72097±6392
|
51945±4690
|
1321.7±54.2
|
6.33±0.06
|
10.95±0.295
|
Bailey
|
n=16
|
92082±19018
|
68099±13872
|
1244.6±47.5
|
6.21±0.04
|
10.86±0.197
|
Randall
|
n=15
|
80300±6162
|
57487±4504
|
1678.1±49.4
|
6.59±0.09
|
11.91±0.241
|
Dallam
|
n=13
|
67783±5025
|
49508±3698
|
1670.0±54.0
|
6.71±0.08
|
12.37±0.271
|
Habitat
|
|
|
|
|
|
|
Rural
|
n=28
|
74488±4147
|
53783±3007
|
1678.9±35.7
|
6.64±0.06
|
12.12±0.182
|
Urban
|
n=42
|
79094±7735
|
57953±5681
|
1337.8±43.5
|
6.33±0.05
|
11.16±0.186
|
Table 3 Coefficients from the linear regression model of predictor variables affecting fecal microbial alpha diversity (observed richness) of black-tailed prairie dogs
Coefficient
|
Estimate (β)
|
SE
|
t-value
|
p-value
|
Intercept
|
1244.62
|
60.81
|
20.469
|
< 0.001
|
Dallam
|
435.30
|
90.82
|
4.793
|
< 0.001
|
Hockley
|
77.09
|
89.01
|
0.966
|
0.3896
|
Lubbock
|
236.46
|
92.88
|
2.546
|
0.0133
|
Randall
|
433.44
|
87.41
|
4.958
|
< 0.001
|
Intercept
|
1337.88
|
38.50
|
34.746
|
< 0.001
|
Rural
|
347.05
|
60.88
|
5.602
|
< 0.001
|
Fecal microbiota beta diversity
The community membership was summarized as Bray-Curtis dissimilarities, unweighted and weighted UniFrac distances significantly explained by sampling sites (F = 3.49, df = 4, p < 0.01; F = 4.33, df = 4, p < 0.01; F = 3.54, df = 4, p < 0.01, respectively) as well as habitat (F = 3.98, df = 1, p = 0.001; F = 5.74, df = 1, p = 0.001; F = 3.95, df = 1, p = 0.001, respectively) in a multiple predictor PERMANOVA model. This test was followed by post hoc pairwise testing between sites using permutational multivariate anlayses of variance using distance metrics (ADONIS) through which it was found that all comparisons were significant (p < 0.05; Supplementary Information: Table S2). The Mantel test showed that geographic separation of samples was correlated with community dissimilarity matrices (Bray-Curtis dissimilarities: r = 0.044, p = 0.003; unweighted UniFrac: r = 0.065, p = 0.024; weighted UniFrac: r = 0.013, p = 0.029). To evaluate the effect of sampling sites and habitat, principal coordinate analysis (PCoA) based on Bray-Curtis dissimilarities, unweighted and weighted UniFrac distances resulted in 16.27%, 22.24%, and 26.72% of the total variation in the community matrix being summarized on the first two axes (Fig. 6A-6C). Multivariate homogeneity of group dispersions (β-dispersion) using Bray-Curtis, unweighted and weighted UniFrac distances corroborated significant variation in the dispersion of samples from group centroids across sampling sites (ANOVA: F = 5.62, p < 0.01; F = 6.12, p < 0.01; F = 6.27, p < 0.01, respectively; Fig. 7A-7C).
Relative abundance of bacterial taxa
To determine which taxa were contributing to the observed differences across sampling sites, we examined differences in relative abundance at the family level. Differences were present in the relative abundance of the family with the most abundant family in Dallam being Bacteroidaceae (22.87%). Prevotellaceae (21.91%) was significantly higher in all samples in Randall County compared to other sites (Fig. 8). Bailey and Hockley were largely dominated by Rikenellaceae and Porphyromonadaceae (26.63% and 15.84% respectively), while Lubbock was mostly abundant with Bacteroidaceae (16.28%). Using a hierarchical clustering method based on the relative abundance of genera could fairly separate the samples from each sampling site (Fig. 9). Two distinct bacterial groups were differentiated by their divergent abundance patterns in the samples. The highly abundant group in the samples including the genera Alistipes, Bacteroides, Prevotella, Roseburia, Parabacteroides, Anaerovorax, Blautia, Eubacterium, and Clostridium. On the other hand, the group displayed the contrasting patterns constituting Flavobacterium, Barnesiella, Paraprevotella, Lactobacillus, Alkaliphilus, Acetivibrio, Catabacter, Fusicatenibacter, Butyrivibrio, Microbacter, and Robinsoniella. The latter group, however, consisted of some genera (Alkaliphilus, Fusicatenibacter, Robinsoniella) causing infection in humans [82-84]. Moreover, Bailey showed higher contributions of Alistipes (18.37%) whereas Parabacteroides (9.15%) mostly prevailed in Hockley. Both Dallam and Lubbock had a greater abundance of Bacteroides (27.74%, 22.21%, respectively). Prairie dogs from Randall County had higher relative abundance of Prevotella (9.38%).
Differential abundance testing
Species-level comparisons were accounted for by differential abundance analyses (ANCOMBC) with 71 species observed to differ across five experimental groups (Fig. 10). To provide a robust understanding of the microbiome differences between prairie dogs living in rural areas and those living in urban areas, we also used ANCOMBC. This analysis revealed that 43 species were differentially abundant between rural and urban habitats. A greater percentage of Clostridium sp., Cohaesibacter haloalkalitolerans, Prevotella sp., and Thermophagus xiamenensis were found in rural areas (Fig. 11). In contrast, the fecal microbiota of urban areas was characterized by higher proportion of Alistipes sp., Anaerovorax odorimutans, Aureibacter tunicatorum, Bacteroides sp., Butyrivibrio fibrisolvens, Christensenella minuta, Clostridium sp., Desulfosporosinus sp., Eubacterium brachy, Fusibacter sp., Kiloniella sp., Mogibacterium pumilum, Peptostreptococcus russelli, Ponticoccus litoralis, Rhizobium sp., Sphingobium sp., Sporobacter termitidis, and Vallitalea pronyensis.
Furthermore, to explore the bacterial species contributing to the important values in the habitat, we used random forest classifier. We found Alistipes shahii had higher mean decrease accuracy (MDA = 3.53; zotu 906, 1401, 1041, 3705, 905, 4211, 745) compared to others and thus had a greater impact on the accuracy of the classification followed by Bacteroides rodentium (MDA = 3.34; zotu 1572, 98, 3354, 4897, 2937, 589, 1351, 4298, 3210, 3787) and Prevotella shahii (MDA = 2.83; zotu 386, 2852, 1893, 1847, 2777, 2113) (Fig. 12).
Relationship among environmental variables, fecal microbial alpha, and beta diversity
To investigate whether environmental factors were associated with the changes in black-tailed prairie dog fecal microbial diversity, we used linear and linear mixed effect models based on AICc. We observed a significant positive association between Shannon diversity index and cumulative precipitation (Pearson correlation: r = 0.42, p = 0.0003; Fig. 13A), but negative correlation with average maximum (r = -0.47, p < 0.001) and average minimum temperature (r = -0.35, p = 0.0033) (Fig. 13B and 13C). No significant correlation was observed between elevation and Shannon index (r = 0.19, p = 0.11). First principal component of PCoA (based on unweighted UniFrac distance) that explained 13.69% variation of the beta diversity was positively correlated with cumulative precipitation (r = 0.69, p < 0.001) and elevation (r = 0.31, p = 0.01) (Fig. 13D and 13E). In contrast, significant negative correlation was observed with both average maximum and average minimum temperature (r = -0.77, p < 0.001; r = -0.59, p < 0.001, respectively, Fig. 13F and 13G). In both alpha and beta diversity models, average maximum temperature was the strongest predictor considered the top ranked model (F(1,68) = 18.97, R2 = 0.22, p < 0.001; F(1,68) = 54.71, R2 = 0.39, p < 0.001, respectively; Table 4).
Table 4 Alpha (Shannon diversity index) and beta diversity (PCoA 1) models based on AICc model selection
Alpha diversity models
|
K
|
AICc
|
ΔAICc
|
Wi
|
Cumulative Wi
|
LL
|
Average maximum temperature
|
3
|
47.50
|
0.00
|
0.39
|
0.39
|
-20.57
|
Average maximum temperature + Average minimum temperature
|
4
|
48.03
|
0.52
|
0.30
|
0.69
|
-19.71
|
Cumulative precipitation + Average maximum temperature
|
4
|
49.61
|
2.10
|
0.14
|
0.82
|
-20.50
|
Cumulative precipitation + Average maximum temperature + Average minimum temperature
|
5
|
50.33
|
2.82
|
0.09
|
0.92
|
-19.70
|
Cumulative precipitation
|
3
|
51.36
|
3.86
|
0.06
|
0.97
|
-22.50
|
Cumulative precipitation + Average minimum temperature
|
4
|
53.41
|
5.90
|
0.02
|
0.99
|
-22.40
|
Average minimum temperature
|
4
|
55.78
|
8.28
|
0.01
|
1.00
|
-24.71
|
Beta diversity models
|
|
|
|
|
|
|
Average maximum temperature
|
3
|
-158.90
|
0.00
|
0.37
|
0.37
|
83.26
|
Average maximum temperature + Average minimum temperature
|
4
|
-157.32
|
0.00
|
0.06
|
0.91
|
83.34
|
Elevation + Average maximum temperature + Average minimum temperature
|
5
|
-156.82
|
1.08
|
0.22
|
0.59
|
83.88
|
Cumulative precipitation + Elevation + Average maximum temperature + Average minimum temperature
|
6
|
-155.94
|
1.96
|
0.14
|
0.73
|
84.63
|
Cumulative precipitation + Average maximum temperature + Average minimum temperature
|
5
|
-155.59
|
2.31
|
0.12
|
0.84
|
83.26
|
Cumulative precipitation + Elevation + Average maximum temperature
|
5
|
-153.65
|
4.25
|
0.04
|
0.95
|
82.29
|
Cumulative precipitation + Average maximum temperature
|
4
|
-152.78
|
5.12
|
0.03
|
0.98
|
80.70
|
Elevation + Average maximum temperature
|
4
|
-152.24
|
5.66
|
0.02
|
1.00
|
80.43
|
Elevation + Cumulative precipitation
|
4
|
-137.35
|
20.55
|
0.00
|
1.00
|
72.98
|
Cumulative precipitation
|
3
|
-136.74
|
21.16
|
0.00
|
1.00
|
71.55
|
Cumulative precipitation + Elevation + Average minimum temperature
|
5
|
-136.03
|
21.87
|
0.00
|
1.00
|
73.48
|
Cumulative precipitation + Average minimum temperature
|
4
|
-134.67
|
23.23
|
0.00
|
1.00
|
71.64
|
Average minimum temperature
|
3
|
-122.29
|
35.61
|
0.00
|
1.00
|
64.33
|
Elevation + Average minimum temperature
|
4
|
-120.36
|
37.54
|
0.00
|
1.00
|
64.49
|
Elevation
|
3
|
-99.22
|
58.68
|
0.00
|
1.00
|
52.79
|