Subspecific variation in gut microbiota of North American bison in a sympatric setting reveals differentially abundant taxa

Abstract

Gut microbiomes play critical roles in host-environment interactions, reflecting habitat and foraging niches. North American bison (Bison bison) subspecies—plains bison (B. bison bison) and wood bison (B. bison athabascae)—exhibit limited genetic variation from historic population bottleneck events, potentially undermining their evolutionary potential. Understanding variation in gut microbiota composition between subspecies may shed light on genetic, phenotypic, and ecological divergence relevant to their adaptive capacities. Using 16S rRNA metabarcoding of fecal samples, we characterized the gut microbiota of both subspecies in the sympatric environment of Elk Island National Park, providing insight into potential phylogenetic gut microbiome divergence. Like other ruminants, the gut microbial community of both subspecies consists primarily of the bacterial phyla Firmicutes and Bacteroidetes. Subspecific classification explained no significant differences in alpha diversity (p > 0.05) in the overall dataset, but has a potentially significant effect on beta diversity (p < 0.05, R² = 0.04). Gut microbiota divergence between subspecies may be driven by differential abundance of specific taxa and associated functional pathways, likely influenced by dietary preferences, ancestral phenotypes, and historical ranges. Our findings support further investigation into diet-microbiome relationships between subspecies in sympatric environments and metagenomic approaches to explore functional differences in the gut microbiome.

Supplementary Information: The online version contains supplementary material available at 10.1186/s42523-025-00451-7.

Keywords: Gut microbiome; North American bison; Phylogenetic divergence; Subspecific variation.

Fig. 1

Map of Elk Island National Park (EINP) with collection areas indicated. Adapted from “Elk Island National Park Site Map” [57]. Wood bison (n = 100) and plains bison (n = 100) fecal samples were collected between October 21–23, 2023

Fig. 2

Relative abundance bar plot of taxa classified by QIIME2 at the family level for plains and wood bison for the top five most abundant phyla. Up to the top 5 taxonomic families within each of the top 5 phyla are represented here, with all others grouped into the category “Other”. Each bar represents the gut microbial profile of an individual bison, with colour corresponding to microbial phylum and shade corresponding to families within their respective phylum

Fig. 3

Effect size of statistically significant, differentially abundant families within plains and wood bison gut microbial community profiles, as determined by ALDEx2 analysis. A negative effect size indicates an increased abundance in plains bison (yellow), whereas a positive effect size indicates an increased abundance in wood bison (blue). Significance was determined within ALDEx2 by a Benjamini-Hochberg corrected, Wilcoxon rank-sum test p-value of less than 0.05

Fig. 4

PICRUSt2 predicted KEGG pathways significantly different (p < 0.05) between plains and wood bison. Differential abundance was determined via ALDEx2, with significance tested by a Benjamini-Hochberg corrected, Wilcoxon rank-sum test. Negative log2fold changes indicate enrichment in the wood bison gut microbial community, whereas positive log2fold changes indicate enrichment in the plains bison gut microbial community

Fig. 5

Alpha diversity assessment of the gut microbial community of plains and wood bison in EINP. (A) Observed ASV richness for plains and wood bison, (B) Shannon Diversity Index for plains and wood bison, (C) Simpson’s Diversity Index for plains and wood bison, (D) Faith’s PD for plains and wood bison, (E) Observed ASV richness for female and male bison of both subspecies, (F) Shannon Diversity Index for female and male bison of both subspecies, (G) Simpson’s Diversity Index for female and male bison of both subspecies, (H) Faith’s PD for female and male bison of both subspecies. Boxplots A-D were generated using the complete dataset of males, females, and individuals of unknown sex for both subspecies. Boxplots E-H were generated using a subset of the data, with individuals of unknown sex removed from the analysis; boxplots have been visually split to highlight differences between subspecies. Asterisks above boxplots indicate statistically significant differences (p < 0.05) as determined by the Wilcoxon rank-sum test. Lines connecting boxplots denote the specific comparisons where significance was observed

Fig. 6

Beta diversity assessment of the gut microbial communities of plains and wood bison in EINP. (A) Bray-Curtis Dissimilarity comparing plains and wood bison, (B) Bray-Curtis Dissimilarity comparing male and female individuals of both subspecies, visually split to highlight differences between subspecies. Ellipses represent 95% confidence intervals. Significance was tested via PERMANOVA, with assumptions checked via PERMDISP. Subspecies was found to have a significant effect on beta diversity (R² = 0.04, p < 0.05), however the PERMDISP was also significant (F = 8.35, p 0.05); perceived differences between subspecies could be due to within subspecies dispersion as opposed to compositional differences in the gut microbial community. The PERMDISP results were not significant between sexes of either subspecies (p > 0.05); sex differences differed significantly between male and female plains bison (R² = 0.03, p < 0.05), but not between male and female wood bison (R² = 0.002, p > 0.05)