Glucocorticoids coordinate changes in gut microbiome composition in wild North American red squirrels

The gut microbiome impacts host health and fitness, in part through the diversification of gut metabolic function and pathogen protection. Elevations in glucocorticoids (GCs) appear to reduce gut microbiome diversity in experimental studies, suggesting that a loss of microbial diversity may be a negative consequence of increased GCs. However, given that ecological factors like food availability and population density may independently influence both GCs and microbial diversity, understanding how these factors structure the GC-microbiome relationship is crucial to interpreting its significance in wild populations. Here, we used an ecological framework to investigate the relationship between GCs and gut microbiome diversity in wild North American red squirrels (Tamiasciurus hudsonicus). As expected, higher GCs predicted lower gut microbiome diversity and an increase in metabolic taxa. Surprisingly, but in line with prior empirical studies on wild animals, gastrointestinal pathogens decreased as GCs increased. Both dietary heterogeneity and an upcoming food pulse exhibited direct effects on gut microbiome diversity, whereas conspecific density and reproductive activity impacted diversity indirectly via changes in host GCs. Our results provide evidence of a gut–brain axis in wild red squirrels and highlight the importance of situating the GC-gut microbiome relationship within an ecological framework.

. Conceptual model demonstrating how ecological factors structure the relationship between glucocorticoids and the gut microbiome. (A) Environmental covariance results in direct effects on both variables and diminishes a detectable effect of GCs on gut microbiome diversity. (B) Ecology does not influence either variable and a direct effect of GCs on gut microbiome diversity is preserved. (C) Ecological factors influence gut microbiome diversity indirectly via host GCs. Note that the three scenarios are not mutually exclusive, such that a combination of direct (A) and indirect (B) effects may result in the appearance or absence of a relationship between GCs and gut microbiome diversity. GCs if conspecific density also increases due to that elevation in food availability 28 . If elevated density results in more frequent social interactions, it may enhance microbiome diversity directly via increased microbial transmission among conspecifics 29 . Such environmental covariance may drive the absence of a relationship between GCs and microbiome diversity in unmanipulated populations ( Fig. 1) [30][31][32] , necessitating an integrative approach to determine how GCs impact the gut microbiome in wild animals.
In this study, we test the hypothesis that ecological factors structure the effects of GCs on gut microbiome alpha diversity in wild North American red squirrels (Tamiasciurus hudsonicus) living in the Yukon, Canada. Red squirrels are highly territorial animals that experience dramatic shifts in food availability and population density as a result of fluctuations in their preferred food source, seeds from white spruce trees (Picea glauca) 33 . Squirrels incorporate other food sources into their diet when seasonally available (e.g., fungi, bark, leaves, flowers) 34 , resulting in changes in dietary diversity that may directly impact gut microbiome diversity. However, spruce seeds comprise the majority of their diet 34 despite their episodic availability. Masting events occur every 4-6 years in white spruce, resulting in the production of a superabundance of cones containing seeds that become available in the autumn. By contrast, few to no cones are available in non-mast years 33,35 . In anticipation of an upcoming spruce mast, squirrels exhibit an extended breeding season and concomitant behavioral changes: territoriality breaks down and conspecific interactions increase due to increased breeding frequency and infanticidal behavior 36,37 . An upcoming spruce mast may thus exert direct positive effects on gut microbiome diversity via more frequent social interactions, which leads to greater horizontal microbial transmission 29 .
Squirrel densities also fluctuate in parallel with food pulses, with densities at their lowest in the months prior to a mast and highest in the spring following a mast 28,38 . Although sociality is expected to increase gut microbiome diversity in group-living animals 29 , this effect may not occur in territorial species 39 . For example, elevated conspecific densities result in increased frequency of long-range territorial vocalizations emitted by red squirrels in our study population 40 , which can in turn reduce interaction frequency by deterring territorial intrusions 39,41 . Indeed, conspecific interactions in squirrels do not appear to vary with density, and the number of territorial intruders has been both negatively correlated 40 and unrelated 41 to density. Given that both actual and perceived increases in density cause GC elevations in red squirrels 28,42 , density may have indirect rather than direct effects on microbiome diversity due to the psychosocial stress of anticipating greater competition 28 .
We predicted to find an overall negative association between host GCs and gut microbiome alpha diversity, along with an increase in pathogenic taxa and taxa involved in host metabolism with increasing GCs. Similar to prior studies in wild animals (Table 1), we focused here on the unidirectional effects of GCs on microbial diversity to understand gut-brain axis function in the absence of an experimental manipulation, and because previously documented sensitivity of GCs to intrinsic and extrinsic factors in our population 28 suggests it may mediate downstream effects on the microbiome. We used a multivariate structural equation modeling approach 43,44 to integrate GCs, the gut microbiome, and exogenous ecological and host biological variables into a single causal network 45 . We tested a set of a priori hypothesized relationships related to the direct and indirect effects of dietary heterogeneity, an upcoming spruce mast, and conspecific density on both GCs and gut microbiome diversity (Fig. S1). We expected to find that dietary heterogeneity and an upcoming spruce mast would have direct positive effects on gut microbiome diversity. Conversely, we predicted that density would have an indirect negative effect on diversity by way of GC elevations. We additionally included biological factors (reproductive activity, sex, age) in our analysis, given their potential effects on GCs and microbiome composition [46][47][48] . We predicted that reproductive activity had positive direct effects on both GCs and gut microbiome diversity, as a result of increased energetic demands 49,50 and conspecific interactions, respectively. We also expected that older age would predict lower microbiome diversity, and that males would exhibit greater microbiome diversity due to travel across territories for multiple mating in the breeding season 51 .

Results
Gut microbiome diversity is negatively associated with glucocorticoids. Both gut microbiome alpha diversity and GCs were highly variable across seasons in each of our sampling years, with gut microbial diversity reaching its maxima during the summer months of July and August ( Fig. 2A), coinciding with increased dietary diversity 34,52 . GCs were highest in early spring (March), with the exception of the mast year of 2010 in which GCs steadily increased across the first part of the year ( Fig. 2A). Consistent with our predictions and in line with prior studies in which GCs were experimentally manipulated (Table 1), GCs were negatively associated with gut microbiome alpha diversity. Individuals with greater GC concentrations exhibited relatively lower taxonomic diversity (i.e., species richness, Chao1: β = − 75.05 ± 25.91, t = − 2.90, P = 0.004; Fig. 2B). Greater GCs were also associated with lower Shannon Indices, a composite measure of species richness and evenness (Shannon: β = − 77.64 ± 36.00, t = − 2.16, P = 0.03, Fig. 2C), as well as lower phylogenetic diversity in the gut microbial community (Faith's PD: β = − 66.97 ± 26.71, t = − 2.51, P = 0.01, Fig. 2D). The negative relationship between GCs and gut microbiome alpha diversity was robust to individual variation in GC production, with higher individuallyaveraged GCs similarly predicting lower species richness (β ± SE: − 0.027 ± 0.011, t = − 2.45, P = 0.02).
The gut-brain axis mediates ecological effects on the gut microbiota. To determine how ecological and host factors contributed to the effects of GCs on gut microbiome alpha diversity, we fit a structural equation model (SEM) based on a set of a priori hypothesized pathways (Fig. S1). The SEM was constructed to www.nature.com/scientificreports/ test the relative direct and indirect effects of three ecological factors (dietary heterogeneity, an upcoming spruce mast, and conspecific density) and three host factors (reproductive activity, age, and sex) on GCs and gut microbiome diversity. Overall, the SEM revealed direct and indirect pathways by which ecological and host factors exert cascading effects on microbial diversity ( Figure 4). In line with our predictions, dietary heterogeneity and the presence of an upcoming spruce mast both exhibited direct effects on gut microbiome diversity, such that a more heterogeneous diet (standardized β = 0.26, P < 0.001) and an upcoming spruce masting event (standardized β = 0.71, P < 0.001) led to greater microbial diversity. Tests of directed separation revealed that there was no effect of either dietary heterogeneity or an upcoming spruce mast on host GCs. By contrast, conspecific density and reproductive activity indirectly, but not directly, affected gut microbiome diversity via changes to host GCs. Higher squirrel densities (standardized β = 0.26, P = 0.006) and reproductive activity (standardized β = 0.47, P < 0.001) predicted greater GC concentrations, which in turn reduced microbial diversity (standardized β = − 0.18, P = 0.01). There was no effect of age or sex on gut microbiome diversity, and tests of directed separation similarly found no effect of age or sex on host GCs.

Figure 3.
Gut microbial taxa that shift with increasing host glucocorticoids. Barplots depict bacterial families (A) and genera (B) whose relative abundance was significantly (Benjamini-Hochberg adjusted P < 0.05) predicted by changes in host glucocorticoid concentrations. Bold taxa exhibited a mean relative abundance > 0.01%. Effects of GCs reflect model estimates generated by negative binomial mixed models testing the effect of GCs on the relative abundance of each bacterial taxa, controlling for collection date, food supplementation status, and individual ID. Black bars represent a decrease in relative abundance with increasing GCs; grey bars represent an increase with increasing GCs. Taxa depicted at the bottom of Panel B (Kribbella, Propionibacterium) exhibited model estimates ~ 10× larger than the rest of the taxa and are therefore separated from the main plot for visualization purposes.

Discussion
Determining if ecological factors structure the relationship between glucocorticoids (GCs) and gut microbiome alpha diversity is crucial to interpreting the adaptive value of the GC-microbiome connection in wild animals.
Here, we demonstrate that the link between GCs and gut microbiome alpha diversity is robust to ecological factors that directly influence the gut microbiome in wild populations. Our findings additionally suggest that GCs can integrate changes in ecology and host biology to induce plasticity in the gut microbiome. Gut microbiome alpha diversity varied across seasons and peaked in summer, coinciding with the period of greatest dietary heterogeneity 34,52 (Fig. 2A). GCs exhibited greater variability overall: GCs reached their peak in the early spring in non-mast years (2008, 2009), but in summer of the mast year (2010) ( Fig. 2A). Despite this variability, there was a significant negative association between GCs and gut microbiome diversity such that individuals with greater GCs exhibited lower alpha diversity across three separate measures (species richness, evenness, and phylogenetic diversity) (Fig. 2B-D). Though prior studies in unmanipulated populations did not detect a relationship between GCs and gut microbiome alpha diversity (Table 1), but see 55 for oral microbiome), our results align with experimental studies in which GCs were manipulated. This suggests that gut-brain axis communication, particularly between GCs and gut microbiome composition, is under strong selection in our population.
As GCs increased, the taxa that decreased in the gut microbiome were overwhelmingly rare taxa that contributed < 0.01% in relative abundance to the overall community (e.g., Odoribacteriaceae, Sporichthyaceae) (Fig. 3). This finding aligns with expectations about the effects of disturbances on gut microbiome composition from an ecological perspective 56 . Resilience to microbial disturbances is greater among abundant taxa than nonabundant taxa 57 , likely due to divergent patterns of colonization and succession. Abundant taxa retain their position in microbial communities via selective filtering and by occupying core niches 56 . By contrast, rare taxa are incorporated into the microbiome largely via stochastic processes, though they contribute significantly to community alpha diversity measures 58 . A reduction in rare taxa with increasing GCs in our population may therefore indicate that the effects of elevated GCs on rare bacteria in the gut mimic those expected by broader microbial disturbances (e.g., antibiotics, infection) 59 .
A decrease in the relative abundance of rare bacteria may also serve to reorganize host metabolic priorities through replacement by core taxa that can better support changes in host energetic demands. Overall, we found that increases in host GCs were accompanied by increases in bacteria involved in host metabolism (Fig. 3). Both Oscillospira, which correlates with the consumption of spruce buds in the late spring 52 , and Coriobacteriaceae, a common rodent gut microbe involved in energy metabolism 60 , increased in relative abundance with increasing GCs. In experimental settings, housing stress caused an increase in Coriobacteriaceae 61 , suggesting that it may similarly contribute to maintaining energy balance in wild rodents facing challenging environmental conditions. We additionally found that individuals with elevated GCs exhibited greater relative abundances of Ruminococcaceae, a bacterial family of cellulolytic and fibrolytic bacteria involved in dietary acclimation in wild animals 20,62 . Together, our results suggest that in our study population, GCs may coordinate adaptive shifts in gut microbiome composition in response to increased energetic demands, seasonal changes in diet, or both.
Resistance to pathogens has been proposed as one of the major evolutionary advantages conferred by host microbial communities 63 and butyrates, compounds produced via fermentation by microbiota in the large intestine 64 , are particularly critical to preventing intestinal pathogen invasion 65 . We found that the butyrateproducing bacteria Butyriciococcus (family Ruminococcaceae) and Clostridium were elevated in the red squirrel www.nature.com/scientificreports/ gut microbiome as GCs increased (Fig. 3). In line with this finding, elevated GCs were associated with lower relative abundances of two potentially pathogenic genera: Salmonella and Yersinia. Salmonella are rod-shaped, Gram-negative bacteria that can cause gastroenteritis in both rodents and humans upon infection 66 . Similarly, Yersinia (including Y. pestis, Y. pseudotuberculosis, and Y. enterocolitica) are commonly harbored in the gut microbiota of wild rodents and can lead to enteric and systemic disease 54,67 , though we note there is currently no evidence of Yersinia disease in our population. While a loss of these taxa with increasing GCs contradicts theoretical expectations of pathogen susceptibility as microbial diversity decreases, our findings align with a prior study in free-living birds in which elevated GCs similarly reduced the relative abundance of intestinal pathogens 24 , and in piglets in which GCs reduced Salmonella, specifically 68 . These data suggest that elevations in GCs may confer short-term protection against intestinal pathogens, potentially through transient increases in gut immune function or butyrate production. Confirmation of the pathogenicity of these taxa, however, requires high-resolution, strain-level genomic data outside the scope of the present study.
To disentangle the effects of ecological and host factors on the relationship between GCs and the gut microbiome, we constructed a structural equation model (SEM) based on a set of causal a priori hypothesized pathways (Fig. S1) 43,69 . As predicted, an upcoming spruce mast had a direct positive effect on gut microbiome alpha diversity, and this path was the strongest path in the SEM (Fig. 4). An increase in gut microbiome diversity in the mast year compared to the non-mast years aligns with our expectation of territorial breakdown and increased social interactions due to an extended breeding season in the months leading up to a masting event 36,37 . Given that the positive link between sociality and gut microbiome diversity is well-supported at least in some taxa 8,29,70 , squirrels may exhibit increased gut microbiome diversity due to greater horizontal transmission of microbes between conspecifics as social interactions become more frequent.
Similar to the effects of an upcoming spruce mast, dietary heterogeneity had a direct positive effect on microbial diversity, though this effect was approximately 2.5× weaker than that of the upcoming mast (standardized β = 0.26). Gut microbiome alpha diversity was greatest in the months in which the food available to red squirrels was most heterogeneous (e.g., fungi, buds, and seeds) 34 . Indeed, a varied diet is expected to increase microbiome diversity through greater substrate selection for diverse ecological niches 71 . This effect of dietary heterogeneity on gut microbiome diversity, coupled with prior work on the relationship between diet and gut microbiome composition in this population 52 , indicates that the red squirrel gut microbiome responds rapidly to shifts in food availability.
By contrast, conspecific density indirectly, but not directly, impacted gut microbiome diversity by way of GCs, lending support to previous findings that the frequency of social interactions (and thus horizontal transmission) is not related to squirrel density in this population 40 . In line with our expectations, elevated conspecific densities predicted increased host production of GCs 28 , which in turn predicted reduced gut microbiome diversity (Fig. 4). Red squirrels are highly sensitive to changes in density, and signals of both actual and perceived elevated density lead to GC increases independent of other ecological factors that covary with it (e.g., food) 28 . That elevated density reduced gut microbiome diversity via increasing GCs aligns with our understanding of the regulation of the gut-brain axis by the social environment, stress, and psychological state in laboratory rodents 72 . In the wild, vocalizations can buffer individuals from physical interactions with conspecifics even in times of high densities in highly territorial animals like red squirrels. Our results suggest that the indirect effects of increased density on gut microbiome diversity likely reflect the psychosocial stress of increased competition, demonstrating a novel link between social stress and the gut-brain axis in a wild mammal.
Reproductive state was the only host factor to exhibit an effect on gut microbiome diversity, and the effect was indirect and the strongest path to GCs in the model. As predicted, being reproductively active (e.g., scrotal for males, and breeding, gestating, or lactating for females) predicted greater GCs, and in turn a reduction in microbial diversity. In both males and females, reproduction increases host metabolic demands broadly 50 and, in females in our population, GCs specifically 49 . A consequent reduction in microbiome diversity may therefore better support host energy balance by increasing the relative abundance of core microbiota at the expense of rare taxa that contribute less to the metabolic functions of the community. Contrary to our predictions, we found no effect of sex or age on microbial diversity, and the SEM did not identify an effect of either on GCs via tests of directed separation 69 . Studies in humans and other mammals have found mixed effects of age on both GCs and microbial diversity 48,73,74 . Sex effects on microbiome composition related to hormones and behavior have been documented in experimental rodent models 47,75 , but studies in wild populations have not typically found sex differences in gut microbiome composition 76 .
A number of potential pathways may explain how host GCs impact gut microbiota. First, GCs can alter lipid metabolism, leading to lipid accumulation in the gut 77 . Increased lipid metabolism reduces taxonomic diversity in the gut microbiome of laboratory rodents as some bacteria exhibit sensitivity to lipid accumulation 78 . Second, elevated GCs can disrupt host circadian rhythms in laboratory settings 78 , reducing gut microbiome alpha diversity 79 . Finally, increased GCs can decrease mucin synthesis, which is integral to the stability of the gut microbiome and largely determines its composition 80 . Impaired mucin synthesis as a result of elevated GCs may therefore disrupt gut microbiome stability by reducing the resilience of the mucosal layer, leading to a reduction in non-core bacteria and an overall loss of community diversity 81 .
Of note is the bidirectionality of the gut-brain axis demonstrated in laboratory rodent studies 82,83 , and thus the potential bidirectionality of the GC-microbiome relationship in wild animals. Here, we focused on the unidirectional effects of GCs on gut microbiome diversity similar to prior studies (Table 1). However, gut microbiota can themselves regulate the hypothalamic-pituitary-adrenal axis, such that shifts in gut microbiome composition may directly modulate host production of GCs 84 and contribute to a feedback loop between the two systems 85 . Statistical constraints inherent to structural equation modeling prevented us from incorporating a bidirectional relationship between GCs and microbial diversity into this study 69 . Nonetheless, the bidirectionality of the gut-brain axis has important implications for the evolution of the relationship between GCs and microbial www.nature.com/scientificreports/ diversity in wild mammals 11 . We encourage future research on wild populations to implement experimental frameworks whenever possible to better characterize the complexity of the GC-gut microbiome relationship.

Methods
Ethics statement. All methods were carried out in accordance with relevant guidelines and regulations. All research methods and protocols were conducted under animal ethics approvals from Michigan State University (AUF#04/08-046-00), University of Guelph (AUP#09R006), and University of Michigan (PRO00005866). All authors complied with the ARRIVE guidelines.
Study population. Subjects for this study were wild North American red squirrels (Tamiasciurus hudsonicus) inhabiting a natural environment in southwest Yukon, Canada (61° N, 138° W). All subjects were continuously monitored as part of the Kluane Red Squirrel Project, a long-term field study that has been conducting a combination of live-trapping, focal behavioral observations, sampling, and experimental manipulations in the area since 1987 38,86 . Individual squirrels are marked with small metal ear tags and unique combinations of colored wire threaded through the ear tags. Individuals are monitored from birth to death in each year of study, roughly from March to October, using live-trapping and behavioral observations 38 . Individuals included in our study lived on one of three grids (Agnes or AG, Kloo or KL, Sulphur or SU). On AG, individuals were supplemented with peanut butter from October to May in each year for a separate experiment focused on experimentally increasing squirrel population density 28,38 . On KL and SU, no food supplementation was provided. All models controlled for food supplementation status given its potential impacts on gut microbiome composition.
Sampling. We collected 227 samples from 88 individuals across three years (2008-2010). When individuals were captured, they were handled such that their unique identity could be determined (by reading ear tags), sexed, and their reproductive condition could be recorded. Fecal samples were collected opportunistically during live-trapping from underneath the traps using forceps. Following capture and handling of squirrels, fresh fecal samples were collected from underneath the traps, kept on wet ice until they could be frozen at − 20 °C within 5 h of collection in the field. Samples contaminated with urine were not collected and all samples were kept at − 20 °C until analysis. We removed one fecal pellet from each sample using sterilized forceps for microbiome sequencing and then used the rest of the sample to measure fecal GC metabolites.
Age. The age of each squirrel was known as individual squirrels were uniquely tagged in their natal nest when they were ~ 25 days of age and age is accordingly recorded at each trapping event 38 .
Reproductive activity. Squirrels were live-trapped regularly and handled using visualization to determine reproduction state. During each trapping event simultaneous with the collection of microbiome and hormone samples, the reproductive state of the individual was determined via abdominal palpation 86 . Males were considered reproductively active if their testes were scrotal, and not reproductively active if testes were abdominal. For females, pregnancy status was assessed via abdominal palpation for fetuses as well as by examining nipple condition. Females were determined to be reproductively active if they were gestating, lactating, or breeding based on nipple condition. We have previously found that females that are reproductively active (pregnant or lactating) have higher fecal GC metabolites than those that were not reproductively active whereas there were no differences in the effects of reproductive activity (presence or absence of scrotal testes) in males 49 .
Dietary heterogeneity. Although red squirrels consume primarily white spruce seeds, they also consume a number of other foods (e.g., spruce bark and needles, willow leaves and buds, fungi, and bearberry flowers) and thus experience varying levels of seasonal dietary heterogeneity 34 . We coded dietary heterogeneity by ranking the availability of these different foods across seasons from greatest (3) to least (1). Samples collected prior to June of each year were ranked as 1, while samples collected in the month of June and late summer (July-August) were ranked as 2 and 3, respectively.
Conspecific density. Densities (expressed as squirrels per hectare) for each grid of the study (KL, SU, AG) were calculated separately for each year (2008,2009,2010) across the dataset using census data. In May of each year, we determined the number of squirrels owning a territory on our study areas using a combination of livetrapping and behavioral observations. Because squirrels are diurnal, regularly exhibit territorial calls, and their territories are visually conspicuous, we were able to completely enumerate all squirrels living in our study areas.

Sequencing and bioinformatics.
Microbiome data used in this study are a subset of previously published data 52 . DNA extraction and sequencing was performed as described in Ren et al. 52 . Briefly, the V1-V3 hypervariable region of the 16S rRNA bacterial gene was amplified using two universal primers: 27F (5′-ARG GTT TGATCMTGG CTC AG-3′) and 534R (5′-TTA CCG CGG CTG CTG GCA C-3′). Samples were barcoded for PCR amplification, pooled, gel purified, and then sequenced on an Illumina MiSeq using 300 bp paired-end sequences. Sequences were then filtered, quality controlled, and reads were successfully merged using QIIME 87 . Chimeras were removed using USEARCH 88 and sequences determined to be non-chimeric by both de novo and reference-based algorithms were retained. Reads were clustered to OTUs using UCLUST 89 with an identity threshold of 97% (genus-level). Mitochondria and chloroplast were removed, and samples were rarefied to 4000 reads per sample. www.nature.com/scientificreports/ Hormone metabolite analysis. The time period from collection in the field to freezing (~ 5 h) did not impact fecal GC metabolites 49 . We measured fecal GC metabolites using previously validated protocols 49,90 . Briefly, samples were lyophilized for 14-16 h, bathed in liquid nitrogen, and pulverized using a mortar and pestle. A subsample (0.05 g) was then extracted using 80% methanol where the samples were vortexed at 1450 RPM for 30 min followed by centrifuging for 15 min at 2500g 49 . The supernatant was then used in an enzymeimmunoassay that employed an antibody that measures GC metabolites with a 5α-3β,11β-diol structure 91 . We have previously validated this assay and shown that the antibody can accurately measure increases in adrenal production of GCs 49 . We have also shown that our measures of fecal GC metabolites are comparable across assays 92 . Using pooled samples that were run repeatedly on different plates (n = 115) in our laboratory show that the estimates of optical density for these pooled samples were highly repeatable (R = 0.851, 95% CI 0.543-0.925). Using a linear mixed-effects model, we partitioned the variance in the optical density recorded for the pooled samples that were run across these different plates and found that most of the variance was due to the sample itself (85.1%) with little of it being explained by intra-assay variation as all samples were run in duplicate (4.9%) or by inter-assay variation (9.9%).

Statistical analysis.
All statistical analyses were conducted in R (v. 3.5.2.) (R. Core Team, 2015). OTU and taxonomy tables were imported into R and merged into a phyloseq object for downstream analyses using the ape 93 and phyloseq 94 packages. All figures were created in R, with the exception of the conceptual model (Fig. 1) and the structural equation model figures (Figs. 4 and Fig. S1), which were created in bioRender (www. biore nder. com).
Alpha diversity. The estimate_richness() function in the phyloseq package was used to calculate the observed richness (Chao1) and Shannon Index of alpha diversity. Faith's Phylogenetic Distances were calculated using the pd() function on the phyloseq object in picante 95 . Linear mixed-effects models were used to assess the relationship between bacterial diversity and fecal glucocorticoid metabolites (GCs), including individual ID as a random intercept, and collection date and food supplementation as fixed effects. GCs concentrations were log-transformed to improve model fit. Shannon Indices were Tukey transformed prior to analysis to achieve residual normality. All models were assessed for multicollinearity among predictor variables by calculating variance inflation factors (VIF < 5).

Differential abundance testing.
To identify the bacterial taxa whose relative abundances were significantly associated with changes in host GCs, we constructed negative binomial mixed models and implemented our analysis using the NBZIMM package 96 . Negative binomial models outperform other traditional differential abundance methods (e.g. DESeq) because they are better equipped to handle the zero-inflation and sparsity common to microbiome count data 96 . Taxa included in differential abundance testing were filtered with a liberal threshold of > 0.001% relative abundance to the overall microbiome community to avoid excluding rare taxa as they contribute substantially to measures of community diversity 58,97 . Models included the read count of each bacterial taxa as the dependent variable, GCs (scaled to zero mean and unit variance) as a fixed effect, controlling for collection date (fixed), food supplementation (fixed), and individual id (random). Taxa whose negative binomial models did not converge due to a high presence of zeroes were modeled instead with zero-inflated negative binomial models using the glmer.zinb() function in the same package (NBZIMM). We controlled the false discovery rate by applying a Benjamini-Hochberg FDR correction to all p-values. Adjusted p-values < 0.05 were considered statistically significant.

Structural equation modelling.
To integrate ecological and host variables into our model framework investigating the relationship between GCs and gut microbiome diversity, we constructed a structural equation model using (SEM) using the piecewiseSEM package 69 . SEM is an effective way to evaluate direct and indirect effects of multiple variables within complex ecological systems 43 . Unlike traditional variance covariance-based SEM, piecewise SEM approaches allow for the inclusion of random effects, the construction of a single causal network from multiple separate models, and the ability to handle small sample sizes and compare models using Akaike information criterion (AIC) 69 .
Using piecewiseSEM, we investigated whether the relationship between GCs and gut microbiome diversity (endogenous variables, i.e., variables of interest) was moderated by host and/or ecological factors (exogenous variables, i.e., variance outside of the model structure). All categorical variables were converted to numeric variables prior to modeling. To build the SEM, we first constructed two component linear mixed-effect models. The first model tested the effects of conspecific density, reproductive activity, dietary heterogeneity, and an upcoming mast on GCs. The second model tested the effects of reproductive activity, dietary heterogeneity, an upcoming mast, age, and GCs on gut microbiome diversity (Chao1 richness). Both component models included sample collection date and food supplementation status as a fixed effect and individual ID as a random effect. However, food supplementation did not affect GCs or gut microbiome diversity in either of the component models (effect on GCs: β ± SE − 4.86 ± 3.08, t = − 1.58, P = 0.12; effect on gut microbiome alpha diversity: β ± SE − 29.14 ± 79.13, t = − 0.37, P = 0.71), and was therefore removed from the SEM to improve model fit (AICc) and refine the standardized beta estimates. The overall fit of the SEM was evaluated using Shipley's test of d-separation Fisher's C statistic and AICc.

Data availability
All sequences, hormone data, and R code related to this manuscript are available at figshare ( www.nature.com/scientificreports/