What this is
- This research examines how environmental quality impacts the gut microbiota of Siberian flying squirrels through hormonal pathways.
- The study focuses on the , particularly how relate to gut microbiota diversity in adults and juveniles.
- Findings reveal that habitat quality affects glucocorticoid levels, which in turn influence microbial diversity, but this relationship varies by age.
Essence
- Environmental quality influences gut microbiota composition in Siberian flying squirrels indirectly through , with age-dependent effects observed. Adults show a positive correlation between and microbial diversity, while juveniles do not.
Key takeaways
- Higher in adult squirrels correlate with increased gut microbial diversity. This suggests that stress responses may enhance microbial richness in adults.
- Juvenile squirrels exhibit lower glucocorticoid levels and lack the same relationship between and gut microbiota diversity, indicating developmental differences in stress responses.
- Habitat quality impacts glucocorticoid levels in adults, with smaller habitat patches leading to higher stress hormones, which subsequently affect gut microbiota composition.
Caveats
- The study's cross-sectional design limits the ability to infer causal relationships between habitat quality, , and gut microbiota.
- Findings may not generalize to other species or ecosystems, as the responses observed are specific to Siberian flying squirrels and their unique habitat requirements.
Definitions
- Gut-brain axis: A bidirectional communication system between the gut microbiota and the brain, influencing physiological and behavioral responses.
- Glucocorticoids: A class of steroid hormones involved in the stress response, affecting metabolism and immune function.
Simplified
Introduction
Anthropogenic change is rapidly reshaping ecosystems with profound effects on organismal fitness and population resilience (Sih et al. 2011). Underlying these effects are the molecular processes that can constrain or catalyse animal responses in both predictable and unpredictable environments (Fuller et al. 2010; Petrullo et al. 2025; Todgham and Stillman 2013). Among vertebrates, the hypothalamic–pituitary–adrenal (HPA) axis is an evolutionarily conserved regulator of the physiological stress response, modulating host responses through a cascade of feedback interactions (Sapolsky et al. 2000). This cascade culminates in the production of glucocorticoid hormones (GCs), including cortisol (the primary GC in humans) and corticosterone (the primary GC in rodents, (Sapolsky et al. 2000)). Acting as an integrator network, the HPA axis can transduce cues of changing conditions and coordinate phenotypic shifts that enhance fitness (Dantzer 2023; Martin et al. 2011). Because the HPA axis is sensitive to a wide array of internal and external inputs, GC concentrations can index changes in host physiology following shifts in diet, predation risk, population density, temperature, and the social environment (Sheriff et al. 2011; Wingfield et al. 1998).
Notably, the HPA axis has both target and non‐target effects on other physiological processes and organ systems. This includes interfacing with the symbiotic community of microorganisms that reside in the gastrointestinal tract (hereafter, gut microbiota), through a secondary 'gut‐brain axis' (Cryan and O'Mahony 2011). Through the gut‐brain axis, neuroendocrine function may drive internal selective pressures that regulate within and among‐host microbial variation (Cryan et al. 2019; Stothart et al. 2019). This variation can be further augmented by external drivers of microbial variation like diet, environment, and host‐to‐host microbial transmission (Baniel et al. 2021; Orkin et al. 2019; Raulo et al. 2021; Ren et al. 2017; Sarkar et al. 2024). In turn, these relationships may influence host fitness: gut microbiota can enhance pathogen protection (Spragge et al. 2023; Stecher and Hardt 2011) and facilitate dietary flexibility (Cani et al. 2019; Kohl et al. 2016), but can also impair immune function and digestive efficiency (Pham and Lawley 2014). Further, gut microbiota can themselves influence HPA axis function (Burokas et al. 2015), and consequently, host behavior (Davidson et al. 2020; Wu et al. 2021). The gut microbiota is thus governed by a series of complex internal feedbacks that span the microbiota and the brain, generating covariation among HPA function, gut microbiota, and the broader environment with many causal linkages (Figure 1).
To date, empirical data has largely supported the existence of a gut‐brain axis in wild animals. For example, synthetic corticosterone administration reduces gut microbial diversity in yellow‐legged gulls (Larus michahellis, (Noguera et al. 2018)). Applications of these and similar findings (e.g., (Blanco et al. 2021)) remain limited, however, in the absence of ecologically‐relevant, naturalistic stressors that may exert context‐dependent effects on the coupling of the neuroendocrine system and gut microbiota (Crespi et al. 2013). Studies in unmanipulated populations are crucial to adding this context (Greyson‐Gaito et al. 2020; Levin et al. 2016; Petrullo et al. 2022). Moreover, theoretical frameworks (e.g., the reactive scope model) suggest that inconspicuous feedbacks between an organism's symbiotic microbiota, neuroendocrine system, and broader ecology impact host stress sensitivity (Houtz et al. 2022) and, in turn, ecological adaptation (Alberdi et al. 2016). The integration of environmental cues by the HPA axis could induce shifts in microbial community composition that promote reconfiguration of gut microbiota to optimally track changing dietary or immune landscapes (Baniel et al. 2021; Orkin et al. 2019). Alternatively, habitat degradation and habitat loss may increase GC production to adaptively mobilise energy stores (Busch and Hayward 2009), but with collateral effects on the microbiota (e.g., constraining taxonomic and/or functional diversity, (Amato et al. 2013; Petrullo et al. 2022; Stothart et al. 2016)).
In this study, we disentangle covariation across the HPA axis, gut microbiota, and environment in a wild mammal. The Siberian flying squirrel (Pteromys volans) is an arboreal specialist species adapted to natural state old‐growth boreal forests and listed as a vulnerable species within the European Union (Santangeli, Hanski, and Mäkelä 2013). Because their dietary niche is relatively narrow, they are particularly susceptible to changes in forest structure. Habitat loss due to commercial logging is a primary threat because it fragments their habitat into variably sized patches of suitable, old‐growth forest separated by unsuitable fields of young commercially grown forest (Santangeli, Wistbacka, et al. 2013; Wistbacka et al. 2018). We test the hypothesis that habitat quality influences gut microbial variation in part through HPA axis activation (i.e., an indirect effect via the gut‐brain axis). We take an integrative approach that includes structural equation modelling to test alternative hypotheses related to direct (through transmission of microbes from environment) versus indirect effects (through environment effects on stress physiology and stress effects on microbiota) without requiring experimental manipulation of a vulnerable species (Eisenhauer et al. 2015; Grace et al. 2010; Laughlin and Grace 2019). We expect that individuals inhabiting poorer quality habitats (smaller patch sizes, lower habitat diversity, and more unsuitable habitat) will exhibit higher GCs, indicating chronic HPA axis activation (Sheriff et al. 2011; but see Santangeli et al. 2019). In line with earlier studies, we predict that elevated GCs will alter gut microbial community ecology through microbial homogenization and a loss of microbial diversity (Noguera et al. 2018; Petrullo et al. 2022; Stothart et al. 2016, 2019).
Causal linkages among the environment, gut microbiota, and hypothalamic–pituitary–adrenal (HPA) axis. The left and right sides depict two different environments, one that is in a natural state and has low levels of disturbance and one that is structured by human action and has high levels of disturbance.
Materials and Methods
Study Species and Field Data Collection
Siberian flying squirrels (Pteromys volans) are arboreal rodents adapted to climax‐state spruce‐dominated boreal forest. Within the European Union, most flying squirrels inhabit Southern Finland, which represents the western edge of their range across Eurasia. As old‐growth forest specialists, flying squirrels are threatened by commercial forestry practices. Here, boreal forests are typically cut before they reach the climax state optimal for flying squirrels. As this is an area of intense commercial forestry, few disconnected patches of natural old‐growth forest remain amidst younger commercial stands. Squirrel populations are thus scattered as a metapopulation across habitat patches of variable quality.
In June 2015, we sampled 78 Siberian flying squirrels (48 adults, 31 juveniles) from nest boxes across 69 sites in Finland. These sites are part of a long‐term study on flying squirrels and are deliberately positioned across a wide range of habitats (e.g., from urban park to natural old‐growth boreal forests varying in their level of forestry‐related disturbance (Santangeli, Wistbacka, et al. 2013; Selonen et al. 2014; Wistbacka et al. 2018)). We went through nest boxes in early June, at a time when flying squirrels are breeding and pups are on average ~6 weeks old. We caught 79 individuals (48 adults and 31 pups). In most cases, nest boxes had either a solitary adult or one adult female and 1–4 pups. We collected data in situ, releasing squirrels to their point of capture < 15 min within capture. Data collection included sexing and aging each individual, recording their body mass and femur length (see (Selonen et al. 2016) for more details), and tagging each individual with a subcutaneous PIT‐tag.
Hair and Faecal Sample Collection
Hair is a minimally invasive sample medium that can index long‐term stress exposure across weeks to months (Gormally and Romero 2020). Hair GCs can thus capture chronic environmental stress more consistently than blood plasma GCs, which are affected by short‐term stressors like handling at the time of collection (Desantis et al. 2016; Sheriff et al. 2011). We clipped a small amount (~10 mg) of hair from the side of the tail of each individual squirrel to measure hair glucocorticoid content as a proxy for long‐term stress exposure of our study subjects (Heimbürge et al. 2019). Hair sampling was done in a systematic and consistent way for all individuals: samples were cut as close as possible to the base of the hair, and all samples were cut in the same exact way, from the same part of the tail, by the same person. This ensured that the GCs measured from the hair were as comparable as possible between individuals, as they would all relate to a similar time span of hair growth during the same time of the year. Hair was stored in a dry place at room temperature before being sent to the Berlin Zoo for further analysis.
Faecal samples were collected directly and opportunistically during animal handling and transferred into clean labelled polypropylene tubes with 1 mL of 98% ethanol. To avoid contamination, tubes and faecal pellets were handled with gloves, which were changed between every individual. Faeces were stored in a fridge for the duration of the field work (2 weeks) and subsequently in a −80°C freezer, until they were shipped frozen to the Earth Microbiome Project at the University of Colorado for sequencing.
Habitat Quality Measures
We used data from yearly habitat surveys processed in a Geographic Information System (GIS) to classify habitats across the study sites into three discrete classes following Wistbacka et al. (2018). Suitable habitat refers to breeding habitat for females, defined as layered spruce‐dominated mixed forests, containing trees of different sizes and ages. Semi‐suitable habitat is defined as monospecific forests (typically pine plantations) where flying squirrels can move but cannot breed because of lack of resources (e.g., food provided by deciduous trees or shelter provided by spruce trees (Wistbacka et al. 2018)). Unsuitable habitat is represented by open areas covered by young sapling stands or clear‐cuts, roads, fields, water, or built‐up areas (i.e., areas unusable for flying squirrels to move, feed, and/or breed). We then calculated three measures of habitat quality: (1) patch size (area in ha), (2) patch diversity (different types of habitats present within 300 m of an individual's nest box) and (3) proportion of suitable area within a 100 m radius from the center of each nesting site, which is where the nest boxes are located. The use of the 100 m radius buffer is based on the known size of the core area used by breeding females and was previously found to be the most relevant radius within which the habitat has the highest influence on the species (Santangeli, Hanski, and Mäkelä 2013; Wistbacka et al. 2018). The extended buffer of 300 m used for patch diversity allows capturing more variance in this metric compared to a 100 m buffer (which would generally have very low diversity and become nearly uninformative). While large, this 300 m buffer still captured the core area used by females during breeding, as well as extended areas and habitats therein in its surroundings. The amount of unsuitable habitat as well as the size of the patch have been found to capture relevant habitat variation within the range typically used by flying squirrels (Hanski et al. 2000; Santangeli, Hanski, and Mäkelä 2013; Wistbacka et al. 2018).
Measurement of Hair Glucocorticoids
Hair samples were washed with 100% methanol for 2 min. Methanol was decanted and the dried samples were powdered in a bead mill. Next, hairs were extracted with 400 μL 90% methanol and aliquots of 20 μL, diluted 1:2 with water and analyzed in the enzyme immunoassay (EIA). Hydroxycorticosterone was quantified using a previously validated enzyme‐linked immunosorbent assay (Voigt et al. 2004). Briefly, the assay used a polyclonal in‐house antibody against cortisol‐21‐hemisuccinate (HS), coupled with bovine serum albumin (BSA) raised in rabbits, used together with the corresponding 21‐HS‐steroid coupled to peroxidase (HRP) as the label for the EIA. Antibody cross‐reactivity with different glucocorticoids was as follows: 4‐pregnen‐11a,17,21‐triol‐3,20‐dione (cortisol) 100%, 4‐pregnen‐11ß,21‐diol‐3,20‐dione (corticosterone) 13.2%, and 4‐Pregnen‐21‐ol‐3,20‐dione (desoxycorticosterone) < 0.1%, respectively.
Standards were prepared in duplicates of 20 μL by 1:1 dilutions in assay buffer ranging from 0.2 – 100 pg/20 μL. All the duplicates, together with 100 μL cortisol‐HRP conjugate in assay buffer (50 mM Na2HPO4/Na2HPO4, 0.15 M NaCl, 0.1% BSA, pH 7.4), were added to a plate, which was then incubated overnight at 4°C and then washed four times with washing solution. One hundred fifty microliters substrate solution was added per well (1.2 mM H2O2, 0.4 mM 3,3′,5, 5′‐tetramethylbenzidine in 10 mM sodium acetate, pH 5.5), and the plates were incubated in the dark for 40 min. Fifty microliters of 4 N H2SO4 were then added to stop the reaction and the samples were measured at 450 nm with a 12‐channel microtiter plate reader (Infinite M 200; Tecan). Hormone concentrations were calculated using the Magellan software (Tecan; see (Finkenwirth et al. 2010)). The sensitivity of the assay was defined as two standard deviations from the signal given by the zero blank and thus was 0.2 pg per well.
MicrobialSequencing and Bioinformatics DNA
DNA extraction, amplification, and sequencing was performed in the Knight Laboratory at the University of Colorado in 2016, following Earth Microbiome Project protocols (see http://www.earthmicrobiome.org/emp‐standard‐protocols/16s/↗). In short, DNA extraction was done using MO Bio PowerSoil kits and PCR was run with primers 515F/806R, targeting the V4 region of the 16S rRNA bacterial gene (Caporaso et al. 2012). Resulting amplicons were sequenced using an Illumina MiSeq. Sequences were clustered into OTUs using the UPARSE pipeline (with de novo clustering using 97% similarity threshold and maximum expected error ≤ 2) and assigned taxonomy using GreenGenes. Further preprocessing was done with the phyloseq package in R (McMurdie and Holmes 2013). Based on sample completeness curves produced with the iNext package (Hsieh et al. 2016), we chose to filter out samples with < 10,000 reads as unrepresentative of the community diversity. Using the same package, we also derived asymptotic alpha‐diversity estimates for each sample. Furthermore, we filtered all OTUs assigned to chloroplasts/Cyanobacteria. Finally, taxon abundances were normalised as proportions within their sample. Microbiome variation was then described through two set of metrics: richness (number of unique taxa) and Shannon diversity (an abundance‐weighted measure of diversity) as measures of alpha diversity, and Jaccard similarity (proportion of shared taxa between two samples) as a measure of beta‐diversity.
Statistical Analyses
Identifying Drivers of Variation in Gut Microbiota
To identify the factors that predicted variation in host gut microbiota, we performed a first marginal permutational analysis of variance (PERMANOVA) on microbiome beta‐diversity (Jaccard similarity values) using the adonis2 function in the vegan package (Oksanen et al. 2007). We first tested the effects of host age, sex, date and family membership on microbiome composition in the full dataset (N = 78). We then ran a second marginal PERMANOVA with family set as a blocking factor, controlling for age and sample date, to assess the effects of the three habitat measures (habitat diversity, patch size, and amount of unsuitable area) on gut microbiota composition. We additionally tested for variation in homogeneity of dispersion as a function of host age using the betadisper() function and a generalised linear mixed‐effects model (family = beta).
Differential Abundance of Microbial Taxa
To identify the taxonomic differences between the gut microbial communities of adults and juveniles, we constructed a series of linear mixed effects models using lme4 (Brooks et al. 2017). Models test how age (fixed factor) predicted the relative abundance of bacterial families (arcsine square root transformed), controlling for host sex, date of sample collection, and nest box/family (random effects). We accounted for multiple comparisons and an inflated false discovery rate for multiple hypothesis testing by adjusting all p values using the Benjamini‐Hochberg correction (Benjamini and Hochberg 2000).
Dyadic Models
To test whether variation in glucocorticoids predicted variation in gut microbiota, we built dyadic Bayesian beta regression models using the brms package in R (R Core team 2015), as previously validated (Raulo et al. 2021, 2024). We used pairwise Jaccard similarity values as our response variable and predicted it with pairwise measures for hair glucocorticoids, with sex similarity (0/1 = same/different sex), age similarity (0/1 = same/different age category), spatial distance (metric distance between nest boxes), and nest sharing (0/1 = same/different nest) added as covariates. All covariates either naturally ranged from 0 – 1 or were scaled to do so to make model estimates for all terms comparable. As Jaccard similarity is a proportional metric bound between 0 and 1, a logit link function was used. We chose to use Jaccard–which captures presence/absence rather than abundance–in this model over another beta‐diversity measure because it better captures potential confounding effects of social microbial transmission (Raulo et al. 2021, 2024). As the modelled values are pairwise measures, and as such not independent of each other, a multi‐membership random effect structure was built into the model to account for the effect of the identity of both individuals in each pairwise comparison. The model was run first for all individual pairs in the data and subsequently separately for adults and pups, omitting age similarity as a covariate.
Effects of Host Traits and Habitat Quality on Gut Microbiota and Glucocorticoids
We used linear mixed‐effects models to examine gut microbial alpha‐diversity (richness, Shannon) as a function of host traits and habitat quality. Shannon indices were Tukey transformed to achieve normality of residuals. A first set of models included host age and sex as fixed covariates, and family membership as a random effect predicting gut microbial alpha‐diversity with hair GCs and measures of habitat quality (patch size, patch diversity, and amount of unsuitable area). For models investigating the effects of habitat quality on gut microbiota, we included sex as a fixed factor, and sampling date and time as random factors given work suggesting diurnal patterns of GC secretion from circulation into hair (Sharpley et al. 2010, 2011).
Testing for Direct and Indirect Effects Using Structural Equation Modelling
We constructed structural equation models (SEMs) to simultaneously test for direct and indirect effects of habitat quality on gut microbiota using piecewiseSEM (Lefcheck 2016). Quantitative structural equation modelling is the gold standard for testing a priori causal hypotheses in complex ecological systems. It allows for the partitioning of direct and indirect (mediator) effects of multiple variables nested within one or more mixed‐effects component models (Laughlin and Grace 2019). The piecewiseSEM package specifically improves upon other more traditional SEM approaches because it accounts for smaller sample sizes through Akaike Information Criterion (AIC) model comparison and permits the inclusion of random effects and thus repeated sampling (Lefcheck 2016). It has also previously been used to elucidate gut‐brain axis connections in wild animals (Petrullo et al. 2022). First, we converted all categorical variables to numeric. To retain power, we built separate SEMs for each age class and tested only the measures of habitat quality and covariates that were found to significantly predict GCs and microbial diversity in our previous models. For each SEM, the first component model included hair GCs (log‐transformed) as the dependent variable, habitat quality (patch size, diversity or amount of unsuitable habitat) and sex as fixed effects, and family membership as a random effect. The second component model included gut microbial alpha‐diversity (richness or Shannon diversity) as the dependent variable, hair GCs as the fixed effect and family membership as a random effect. Component models were compiled into one, and structural equation models were run using the psem() function. Model fit was assessed using AIC(c), and standardised coefficients were extracted to allow for comparison of model estimates across the component models (Lefcheck 2016).
Results
Intrinsic and Extrinsic Drivers of Gut Microbial Variation in Flying Squirrels
The Siberian flying squirrel gut microbiota was composed of bacteria from across 14 phyla, with most microbiota from the phyla Firmicutes, Bacteroidetes, Actinobacteria and Proteobacteria (Figure 2A). Within these phyla, samples contained microbiota from 97 families (Figure S1). Specifically, Lachnospiraceae (mean relative abundance = 39.4% across samples), Ruminococcaceae (21%) and Bacteroidaceae (20.4%) were the most highly abundant families, with all other families contributing < 5% to the total gut microbial community (Figure S1).
A marginal PERMANOVA on Jaccard similarity revealed that an individual's nest box, which reflects kinship/family membership, explained a substantial amount of variation in gut microbiota composition (R2 = 0.47, F = 1.17, p = 0.001). Additionally, host age explained some of the variance (R2 = 0.03, F = 1.89, p = 0.001, Table S1, Figure 2B). There was no effect of host sex on gut microbial composition (Table S1). Compared to adults, juvenile squirrels exhibited significantly greater compositional variability (greater distances to group centroid in Jaccard compositional space; estimate ± SE: 0.14 ± 0.03, z = 4.84, p < 0.001, Figure S2, Table S3), suggesting greater inter‐individual differences in microbial composition during early development compared to later life. At the taxonomic level, juveniles housed significantly higher relative abundances of gut microbiota from the families Bacteroidaceae, Rikenellaceae, and Erysipelotrichaceae (Figure 2C, Table S4). By contrast, adults housed greater relative abundances of Coriobacteriaceae, Syntrophomonadaceae and Eubacteriaceae (Figure 2C, Table S4).
Age effects on microbial community composition. (A) The flying squirrel gut microbiota was dominated by microbiota from four main phyla. The most abundant microbial families (ranked) are shown below each phyla in the legend. (B) Principal coordinates analysis of Jaccard similarity across all faecal samples ( = 78) by age class. (C) Differential relative abundances of bacterial families (> 0.05) as a function of age class. Lollipop plot depicts effect sizes from generalized linear mixed‐effects models testing the differential relative abundance of bacterial families (relative abundance > 0.05%) as a function of age class. N P FDR
Glucocorticoids Predict Variation in Gut Microbial Diversity and Similarity
For individuals from which we had paired data on hair GCs (N = 43), we found that juveniles exhibited substantially lower GCs than adults (linear mixed‐effects model predicting hair GCs, estimate ± SE: −0.79 ± 0.13, t = −6.24, p < 0.0001, Table S5, Figure S3) as well as lower gut microbial alpha‐diversity (linear mixed‐effects model predicting microbial richness (number of OTUs), estimate ± SE: −327.81 ± 36.99, t = −8.86, p < 0.0001; predicting Shannon Index, −40.87 ± 4.76, t = −8.59, p < 0.0001; Table S6, Figure S3). Covariation between GCs and microbial diversity was also age‐dependent: among adults, higher GCs predicted greater gut microbial Shannon diversity (linear mixed‐effects model, estimate ± SE: 0.48 ± 0.02, t = 28.61, p < 0.0001), with a marginal effect on microbial richness (estimate ± SE: 26.82 ± 14.02, t = 1.91, p = 0.07, Table S7, Figure 3A). There was no relationship between GCs and either measure of microbial alpha‐diversity among juveniles (Table S7, Figure 3B).
To further parse the relationship between GCs and gut microbiota, we modelled microbial beta‐diversity (Jaccard similarity) as a function of similarity in GCs, controlling for spatial distance, sex and nest sharing (Raulo et al. 2021, 2024). Adult squirrels with more similar GC concentrations shared a greater proportion of microbial taxa than those with more divergent GC concentrations (dyadic Bayesian generalised linear‐mixed effects model, posterior mean = −0.36, CI = −0.60 to −0.12; Figure 3C, Table S8). Additionally, adults in closer spatial proximity to one another also had more similar gut microbiota (posterior mean = −0.23, CI = −0.36 to −0.11; Figure 3C, Table S8). By contrast, among juveniles, neither spatial distance nor differences in GC concentrations predicted gut microbiome similarity (Table S6). Instead, similarity among juvenile microbial communities was predicted only by nest box, such that juveniles shared a greater proportion of their gut microbial taxa with their nest‐mates (i.e., siblings, posterior mean = 1.47, CI = 0.89–2.07; Figure 3D, Table S9). Sex (same versus different sexes) did not predict microbial similarity in either age class.
Covariation between glucocorticoids and gut microbiota is present in adults but absent in juveniles. (A) Adult GCs were positively associated with microbial alpha‐diversity (species richness). (B) Among juveniles, there was no relationship between GCs and microbial alpha‐diversity. (C) Adult squirrels with less similar hair GC concentrations and larger spatial distances between them had less similar microbiota (i.e., effect on the negative side, credible intervals do not overlap 0). (D) Among juveniles, differences in GCs did not predict gut microbiome similarity, but juveniles residing in the same nests shared significantly more gut microbial taxa than juveniles residing in separate nests. Posterior means (points) and their 95% credible intervals (lines) are plotted from dyadic Bayesian GLMM with pairwise microbiota similarity among squirrels (Jaccard similarity) as the response.
Age‐Dependent Environmental Effects on theAxis, but Not the Gut Microbiota HPA
Habitat quality was associated with variation in hair GCs, but the magnitude and direction varied as a function of age class and habitat measure. Adults inhabiting smaller patches had higher GCs (estimate ± SE: −0.15 ± −0.07, t = −2.16, p = 0.038), with the opposite trend in juveniles (Figure 4A, Table S10). Habitat diversity positively predicted GCs in juveniles (estimate ± SE: 0.34 ± 0.14, t = 2.34, p = 0.026) but not in adults (Figure 4B, Table S10). The amount of suitable habitat had no effect on GCs in adults or juveniles (Figure 4C, Table S10). Moreover, none of the three habitat measures had direct effects on gut microbial diversity (richness or abundance‐weighted Shannon diversity) in either age class (Table S11).
Adults and juveniles exhibit differential hormonal sensitivity to measures of environmental quality. (A) Patch size and (B) habitat diversity explain variation in GCs in an age‐dependent manner, with no effects of (C) amount of unsuitable habitat on GCs in either adults or juveniles. Scatterplots show results from linear mixed‐effects models. Points reflect samples (partial residuals); lines indicate regression lines from generalised linear mixed‐effects models.
Environmental Quality Shapes Gut Microbiota Indirectly via Glucocorticoids
Because habitat quality predicted GCs, and GCs had a positive effect on gut microbial richness among adults, we used structural equation modelling to test the hypothesis that habitat‐induced variation in GCs mediated downstream effects on gut microbial diversity. Among adults, smaller patches increased GCs (standardised β = −0.39, p = 0.015), which in turn increased microbial diversity (standardised β = 0.39, p = 0.046, Figure 5A, Table S12). Among juveniles, there was no relationship between GCs and microbial diversity. However, a more diverse habitat increased GCs (standardised β = 0.85, p = 0.001), but had no subsequent effects on microbial diversity (Table S12, Figure 5B).
Age‐dependent and indirect effects of habitat quality on gut microbiota. Among adults, (A) smaller patch sizes increased richness via elevations in GCs. Among juveniles, (B) elevations in GCs induced by habitat diversity had no subsequent effect on richness in pups. Solid arrows indicate significant paths; dashed arrows indicate non‐significant paths. Covariates not shown for visualisation purposes (see Table for details). Model estimates for significant paths are shown with significance levels (* < 0.05, ** < 0.01, *** < 0.001) and reflect standardised beta weights from mixed‐effects component models controlling for sex, sample date, time and nest box. S12 p p p
Discussion
Environmental change can activate both adaptive and maladaptive physiological responses associated with the secretion of glucocorticoid hormones by the HPA axis (Dantzer et al. 2014; Jaimez et al. 2012; Romero 2004; Sheriff et al. 2011). Here, we show that flying squirrels responded hormonally to variation in habitat size, such that hair GCs were higher among squirrels inhabiting smaller patches. Many arboreal animals, including tree squirrels and flying squirrels, are adapted to mature forest habitat (Selonen and Mäkeläinen 2017). Changes in forest structure that reduce the size of usable habitat can lead to fragmentation and resource scarcity that elevate GCs (small mammals: (Boyle et al. 2021), nonhuman primates: (Carlitz et al. 2016; Rangel‐Negrín et al. 2009)). Smaller patches can also increase proximity to extrinsic sources of stress (e.g., humans, predators) at habitat edges, increasing GC production (e.g., squirrel gliders, (Brearley et al. 2010, 2012)). Elevated GCs among squirrels inhabiting smaller patches may thus reflect psychosocial stress, increased vigilance, and/or anti‐predator behaviours (Creel et al. 2009; Voellmy et al. 2014). Moreover, when habitats shrink, squirrels may mobilize energy stores to cope with the increased metabolic burden associated with travelling greater distances for food or mates (Crespi and Denver 2005). These responses may lead to chronic HPA axis activation, causing increased accumulation of GCs into hair (Russell et al. 2012).
Unlike the HPA axis, we found no evidence for direct effects of any measure of habitat quality on gut microbiota. Ambient temperature and vegetation explain gut microbial variation with longitudinal sampling in similar populations, suggesting our cross‐sectional analyses may have been temporally constrained (Liu et al. 2019). Nonetheless, we found support for indirect, age‐dependent environmental effects that run through the HPA axis. Unlike the reduction in microbial diversity documented in other wild animals (e.g., (Petrullo et al. 2022; Stothart et al. 2016)), adult squirrels with higher GCs housed a richer gut microbiota with greater abundance‐weighted alpha‐diversity. At the highest level, GCs may mechanistically enhance microbial diversity via stress‐related behaviours like social grooming in highly social species like red‐bellied lemurs (Raulo et al. 2018). This is an unlikely explanation for microbial diversification in less social systems like Siberian flying squirrels, though behavioural changes mediated by GCs may still play a role: higher GCs could promote increased foraging or vigilance, which may increase opportunities for social and environmental microbial transmission (Voellmy et al. 2014). At the molecular level, microbial diversification may occur via bacterial sensitivity to GC‐induced changes in the biochemical environment of the gastrointestinal tract (e.g., pH, barrier function), which can 'fight' or 'feed' different microbiota (Karl et al. 2018). The immunomodulatory capacity of GCs may diversify microbiota through selective targeting of gut bacteria by the mucosal immune system, namely via secretory immunoglobulin A (sIgA). A change in microbial targeting by sIgA could reduce dominance of 'core' microbiota and open niche space for rarer and less abundant bacteria, increasing alpha‐diversity (Catanzaro et al. 2019). Such sIgA‐mediated effects may be further amplified seasonally as host diets change and microbial diversification is needed to maintain the metabolic flexibility that facilitates dietary transitions, or as metabolic rates change to meet thermoregulatory demands (Han et al. 2023).
While a diverse microbiota is typically viewed as beneficial for hosts, diversity does not reliably index better host health or fitness (Williams et al. 2024). Diversification resulting from elevated GCs could reflect a loss of functional homogeneity generated by canalised microbial responses to ecological cues (e.g., specialised diets). If host HPA axis activation acts as an internal perturbation, it could cause microbial 'scatter' that increases alpha‐diversity, similar to the Anna Karenina principle (Zaneveld et al. 2017). A diversified microbiota could thus have potentially deleterious effects for hosts if it interrupts the adaptive calibration of microbial communities in response to internal or external cues (Petrullo et al. 2025). Yet, adult squirrels had more similar microbiota (i.e., lower beta‐diversity) if their GC concentrations were also more similar, suggesting that microbial communities do not vary idiosyncratically in response to HPA axis activation. Instead, microbial diversification in the footsteps of HPA axis activation may reflect a coordinated suite of adaptive physiological changes to stress. Nonetheless, the directionality and magnitude of these effects are likely to vary across species, resulting in apparent discrepancies in GC‐induced enhancement or depletion of microbial community diversity (Davidson et al. 2020).
Our earlier work found no relationship between habitat quality and GCs, but we did not test for age dependence in these patterns (Santangeli et al. 2019). Adult flying squirrels interact with more dimensions of their environment more directly than juveniles, which are largely bound to their nests where they are exposed to micro‐local phenomena and sibling competition. We found that juveniles exhibited divergent hormonal sensitivity to environmental conditions relative to adults: those inhabiting more diverse habitats exhibited higher GCs, and this was the only environmental variable predictive of the HPA axis during development. A diverse habitat is typically viewed as an indicator of a healthy ecosystem (Rapport 1989), and this could translate into higher GCs if juveniles are exposed to a greater diversity of predators, competitors, or nest parasites. Overall, that the hormonal response was divergent between adults and juveniles suggests differential sensitivity to environmental stress as a function of life history stage (Crespi et al. 2013). In line with this interpretation, while habitat disturbance increases host glucocorticoids in adult mammals (e.g., reviewed in (Dantzer et al. 2014)), this effect becomes inconsistent when juveniles are included in such analyses (Pérez‐Ortega and Hendry 2023).
Because both gut microbiota and the HPA axis exhibit independent and age‐dependent sensitivity to intrinsic and extrinsic inputs (de la Cuesta‐Zuluaga et al. 2019; Sapolsky and Meaney 1986), it is perhaps unsurprising that their coordinated function as an integrated axis also varied with host age (Scott et al. 2017). Indeed, juveniles housed a substantially different gut microbiota than adults, enriched in bacterial families involved in the processing of high‐fat diets like milk during early life in mammals (e.g., Bacteroidaceae, Rikenellaceae, (Comstock 2009)). Juveniles also exhibited lower GCs than adults. Because HPA axis responses are typically attenuated during early development in mammals (Sapolsky and Meaney 1986), GC increases may not be robust enough to influence microbial composition at the community level in juveniles (i.e., dose‐dependency). Unlike adults, juveniles are also continually exposed to conspecific microbiota in the nest: we found that nest box was the only significant predictor of juvenile pairwise microbial similarity. Social microbial transmission across littermates may eclipse glucocorticoid‐induced differences in gut microbiota across individuals through a swamping effect, particularly if the pathways linking microbiota to the host endocrine system have not yet been sufficiently primed (Burokas et al. 2015). Further, maternal transfer of antibodies like sIgA to juveniles via milk may block or offset microbial change mediated by the HPA axis by augmenting endogenous juvenile sIgA production in the gut (Ding et al. 2021). Finally, GC accumulation in the hair shaft is an aggregate measure of host GC production over an extended period of time. It may thus be too coarse a measure to capture brief but pronounced hormonal responses if these responses are more variable in juveniles than in adults (Gormally and Romero 2020). Ultimately, because developmental decoupling of the HPA axis and gut microbiota during early life may explain, at least in part, age‐related differences in gut microbiota (Heitlinger et al. 2017), we encourage future gut‐brain axis research to consider host life history stage.
As interest in how the gut‐brain axis can inform animal ecology and evolution grows, so do the number of outstanding questions. For example, which dimensions of host reproduction and longevity are impacted by microbial responses to ecologically‐triggered hormonal responses? Structural equation modelling offers one approach for capturing gut‐brain axis relationships in wild, vulnerable populations in which experimental manipulations are unethical or impossible (Lane 2006). However, it cannot yet account for the bidirectionality inherent to the gut‐brain axis (Burokas et al. 2015), as gut microbiota can themselves produce hormones and neurotransmitters that influence central nervous system function (Carabotti et al. 2015). Controlled field‐based experiments and omics approaches remain necessary next steps to isolate the molecular pathways linking gut microbiota to the brain within a broader ecological context (Davidson et al. 2020; Greyson‐Gaito et al. 2020).
Author Contributions
Lauren Petrullo carried out the statistical analyses of the data and led the writing of the manuscript. Aura Raulo conceptualised the hypotheses herein, led the data collection, developed statistical methods for the analysis of the data and helped with writing the manuscript. Andrea Santangeli participated in the data collection, derived the habitat diversity measures from GIS data and provided feedback on the manuscript. Ralf Wistbacka developed the field system, led the data collection and provided feedback on the manuscript. Arild Husby helped develop the analysis methods and provided feedback on the manuscript.
Conflicts of Interest
The authors declare no conflicts of interest.
Gut-Brain Axis