Individual variation in the avian gut microbiota: The influence of host state and environmental heterogeneity

The gut microbiota have important consequences for host biological processes and there is some evidence that they also affect fitness. However, the complex, interactive nature of ecological factors that influence the gut microbiota has scarcely been investigated in natural populations. We sampled the gut microbiota of wild great tits (Parus major) at different life stages allowing us to evaluate how microbiota varied with respect to a diverse range of key ecological factors of two broad types: (1) host state, namely age and sex, and the life history variables, timing of breeding, fecundity and reproductive success; and (2) the environment, including habitat type, the distance of the nest to the woodland edge, and the general nest and woodland site environments. The gut microbiota varied with life history and the environment in many ways that were largely dependent on age. Nestlings were far more sensitive to environmental variation than adults, pointing to a high degree of flexibility at an important time in development. As nestlings developed their microbiota from one to two weeks of life, they retained consistent (i.e., repeatable) among‐individual differences. However these apparent individual differences were driven entirely by the effect of sharing the same nest. Our findings point to important early windows during development in which the gut microbiota are most sensitive to a variety of environmental drivers at multiple scales, and suggest reproductive timing, and hence potentially parental quality or food availability, are linked with the microbiota. Identifying and explicating the various ecological sources that shape an individual's gut bacteria is of vital importance for understanding the gut microbiota's role in animal fitness.


| INTRODUC TI ON
The gut microbiome-the enteric microbial community and their genes-is strongly influenced by the environment (Rothschild et al., 2018) and plays a major role in host ecology and evolution (Dethlefsen et al., 2007;Koskella et al., 2017). Although much progress has been made in understanding the interplay between gut microbiota and host phenotypes, research in this field has been restricted primarily to laboratory-based model organisms (Colston & Jackson, 2016;Davidson, Raulo, & Knowles, 2020). In nature, animal's microbiota are far more variable and diverse than their captive counterparts because they are exposed to a wider variety of conditions (Hird, 2017;McKenzie et al., 2017) and because wild hosts are more genetically, and therefore phenotypically diverse than laboratory- The gut microbiota can directly affect the development of key physiological and biological traits of the host (Cox et al., 2014;Hansen et al., 2012;Sudo et al., 2004). Although most of our understanding of these processes comes from mammalian germ-free laboratory models, disruption to microbial gut colonization at critical developmental periods can affect host immunity, metabolism and physical development in other systems such as bees, birds and frogs, (Knutie et al., 2017;Schwarz et al., 2016;Simon et al., 2016;Warne et al., 2019), thus representing a broader taxonomic scope with more ecologically valid microbiota than traditional model organisms (Colston & Jackson, 2016;Pascoe et al., 2017). Cavity nesting birds provide an ideal system for microbiota studies in a natural setting as nestlings can be sampled and monitored easily throughout the early stages of development. Previously we showed in a wild-bird population that nestling gut microbiota diversity predicts nestling weight (Davidson et al., 2021), itself a key predictor of offspring survival and recruitment in great tits (Parus major) (Both et al., 1999;Monrós et al., 2002) and in many vertebrates generally (Bowen et al., 2015;Ringsby et al., 1998;Stige et al., 2019). It is well documented that the gut microbiota changes in early development across taxa (Derrien et al., 2019;Grond et al., 2017;Kohl et al., 2019;Robertson et al., 2019;Stephens et al., 2016). In humans it is assumed that greater gut microbiota variation in infant populations compared to adult populations reflect a shift towards greater microbiota stability (Lozupone et al., 2012), yet the role that ecology plays in this process remains poorly understood.
Many key drivers of ecological variation have the potential to affect the gut microbiota. Exposure to different environmental pools of microbes are probably important, such as those in the air, soil or diet (Grond et al., 2017;Kartzinel et al., 2019;Liddicoat et al., 2020;Ren et al., 2017). The link with diet is especially likely to explain why gut microbiota can vary with habitat (Drobniak et al., 2021;Goossens et al., 2022;Teyssier, Rouffaer, et al., 2018). Vegetation and diet can also vary dramatically within habitats, for example with respect to edge effects (Chen et al., 1992;Wilkin et al., 2007), although it is unknown whether individual microbial communities are affected by these fine-scale processes. Between-seasonal dietary change has been linked to temporal variation in the gut microbiota (Baniel et al., 2021;Hicks et al., 2018;Maurice et al., 2015;Ren et al., 2017).
Within-season variation in the gut microbiota is also probably important because in many animals reproductive timing is critical to ensure breeding coincides with seasonal peaks in food abundance (Brommer et al., 2002;Grüebler & Naef-Daenzer, 2010;Rodríguez et al., 2016;Rubenstein & Wikelski, 2003). Although never tested in the wild, the timing of reproduction could therefore dictate, not just the parent's own microbiota, but also that of their nestlings' through dietary microbes.
Similarly, indirect effects mediated by the host, for example chronic stress, may have knock-on effects on gut microbiota diversity or composition (Noguera et al., 2018;Stothart et al., 2019).
Exercise intensity (Mailing et al., 2019) and energetics (reviewed in Lindsay et al., 2020) are associated with the microbiota, while stress and the microbiota are bidirectionally linked (Sudo, 2014;Sudo et al., 2004). In adults, brood size is an indicator of parental effort (Pettifor et al., 2001) and alters glucocorticoid concentrations (Bonier et al., 2011). While in nestlings, larger brood sizes lead to increased competition for food among siblings (Smith et al., 1989) and higher glucocorticoid hormone concentrations (Greggor et al., 2017;Smith et al., 1988). Thus brood size could influence the microbiota of the parents and their offspring through energetics or stress, though no link has been found to date (Liukkonen et al., 2022). Additionally, the gut microbiota in adults is mediated by a variety of sex hormones (Mallott et al., 2020), and can lead to sex-specific associations between the gut microbiota and environmental sources of variation (Org et al., 2016).
Thus, the composition of gut microbes are expected to vary across time in response to environmental variation and with respect to a whole variety of individual state variables. Nevertheless, individuality and stability in microbiota can also be important, as observed in studies of human health (Chen et al., 2021;Fassarella et al., 2021). various ecological sources that shape an individual's gut bacteria is of vital importance for understanding the gut microbiota's role in animal fitness.

K E Y W O R D S
early development, environmental variation, gut microbiome, gut microbiota, host state, repeatability Consistent among-individual differences indicate limited flexibility, that is, constraints caused by an underlying causal factor-for example, genetic, epigenetic, or early environment effects (Rajilić-Stojanović et al., 2013;Zoetendal et al., 2001). Typically stability in phenotypic variation among individuals is quantified by taking multiple measures from the same individual over time, and estimating repeatability in a mixed model framework (Nakagawa & Schielzeth, 2010). Mixed models also enable an examination of whether individual differences are independent of other factors that can be included as random or fixed effects. Despite its widespread use in animal personality and cognition literature (Bell et al., 2009;Cauchoix et al., 2018), to our knowledge this approach has only been applied once to identify the sources of individual level variation in the gut microbiota (Risely et al., 2022).
Repeatability analysis represents a novel technique that will allow microbiome researchers a greater understanding of microbiome individuality than has been achieved to date (Mallott, 2022).
We explored how the gut microbiota varied with the environment, development, and intrinsic factors using a generalist passerine bird species breeding in nest boxes across heterogenous and fragmented woodlands. The great tit is a widely used model organism in ecological and behavioural studies (Cole & Quinn, 2012;Drent et al., 2003;O'Shea et al., 2018;Tinbergen & Boerlijst, 1990). Given the gut microbiota is correlated with the fitness of great tits in our population (Davidson, Somers, et al., 2021), understanding how the microbiota varies with their ecology is an important next step in this system. Within a single breeding season, we simultaneously investigated how the gut microbiota were related to host state (age, sex, reproductive effort) and a range of environmental factors (habitat type, distance to forest edge, number of siblings in the nest), and whether these explained consistent between-individual differences in the gut microbiota. We had clear a priori reasons to expect environmental and life history traits to covary with gut microbiota variation, and to do so differentially across life stages because the developing gut microbiota in young animals are generally more sensitive to environmental variation than in the more established adult gut microbiota (Derrien et al., 2019;Grond et al., 2017). However, we had no a priori predictions about the direction of such effects in the context of gut microbiota community metrics (i.e., high vs. low diversity and relative taxonomic abundance), owing to the complexity of host-microbe interactions (Douglas, 2018;Foster et al., 2017;Zaneveld et al., 2017). Finally, we conducted a range of repeatability analyses at the individual, the nest, and the woodland site levels to determine whether individual nestlings differed consistently over a short but critical period of their development, from 8 to 15 days old. We further examined whether any consistency at the individuallevel was driven by intrinsic differences among individuals, or instead by shared nest or woodland effects, or indeed a range of other fixed effects as discussed above. Overall, we hypothesise there may be diverse and complex relationships between an individual, their gut microbiota and the environment. Our analyses aim to resolve these interactions to inform future experimental work and our understanding of how the microbiota mediates ecological effects on host adaptation.

| Field monitoring and microbiota sampling
Birds were sampled from nine small woodland fragments across Co.
Cork, Ireland, five of which were mixed/deciduous and four coniferous woodlands (see O'Shea et al., 2018). We collected 262 faecal samples from 204 great tits from 63 nests (see Table S1 for final sample sizes post bioinformatic processing) for 16S rRNA gene sequencing. Nest boxes were monitored during April-June 2016 to determine lay dates, hatching dates and nestling survival. Typically, all individuals in a nest were sampled except for cases where nestlings did not produce a sample within 10-15 min of being placed in sampling apparatus and risked becoming chilled or weak. Individual nestlings were sampled when they were 8 days old (+/− 1 day), and again, if they survived, when they were 15 days old (D8 and D15 birds respectively), at which point parents were also sampled. Birds were placed individually into sterile holding bags inside a heated holding case and naturally-produced faecal samples were collected.
Urea has the potential to affect downstream sequencing (Khan et al., 1991) and was removed through absorption by coffee filters placed as lining in the sampling bags. The faecal matter was collected within 15-20 min of placing birds in the sampling bag, after which birds were returned to the nest. Faecal sacks were ruptured immediately using a sterile inoculation loop and placed in a microcentrifuge tube containing 500 μL of 100% ethanol. Samples were stored at -20°C within 8 h of collection until DNA extraction. D8 nestlings were weighed, and individually identified by clipping the tip of one of their nails, taking care to avoid the blood vessel. Nestlings were again weighed at D15 and ringed with a unique identifiable metal ring (British Trust for Ornithology). We had "repeat samples" from 41 individual nestlings on D8 and D15. Both parents were trapped on the nest where possible and, if not already ringed, were fitted with a British Trust for Ornithology ring, weighed, and aged as either "immature" (first year breeding) or "mature"(second year/+ breeding) and sexed using plumage indicators.

| Illumina MiSeq sequencing
Full library preparation details are described in Supporting Information of Davidson et al. (2021). Briefly, the V3-V4 variable region of the 16S rRNA gene was amplified from the DNA extracts using the 16S metagenomic sequencing library protocol (Illumina).
The DNA was amplified with primers specific to the V3-V4 region of the 16S rRNA gene which also incorporates the Illumina overhang adaptor. Samples were sequenced on the MiSeq sequencing platform (Clinical Microbiomics, Denmark), using a 2 × 300 cycle kit, following standard Illumina sequencing protocols. Four negative controls were carried through to sequencing, two extraction controls (qubit results after extraction: too low to read) and two PCR controls. These samples were subsequently excluded from any analyses following bioinformatic processing. These negative controls indicated there was a small number of potentially contaminant amplicon sequence variants (ASVs) (<0.01% of ASVs) detected by the decontam package using the frequency method (Davis et al., 2018). Therefore, our sequence data may contain some ASVs from the external environment though this should not have systematically biased our results as samples were randomized across plates.

| Bioinformatics
The DADA2 pipeline was used to process the raw sequencing data (Callahan et al., 2016) in R version 3.5 (R Core Team, 2019). Sequence quality was visually inspected. Sequences were trimmed to remove adapters and lower quality reads (median quality scores below 25-30 threshold) at the extremities of the sequence and filtered to remove sequences with expected errors above 1. Read errors were estimated before dereplication. Forward and reverse reads were merged to construct "contig" sequences, which were used to construct a sequence table of ASVs, which in turn counts the number of times each unique sequence is detected. The previous steps were performed for each run separately. Then the separate sequence tables were merged and chimeras removed using the "consensus" method. Taxonomy was assigned to each ASV by RDP's Naive Bayes Classifier (Wang et al., 2007) against the Silva reference database (version 132) (Quast et al., 2012).
Alpha diversity (both Shannon and Chao1 diversity) was calculated using the "estimate_richness" function from the phyloseq package on the filtered data set.
After removing low read samples, chloroplast sequences and mitochondria sequences, there were 18,890,006 total reads clustered into 54,343 ASVs in 246 samples (see Table S1 for sample breakdown). Reads ranged from 10,220 to 557,336, with a mean of 76,789 reads per sample. For the analyses presented here the reads were not rarefied prior to alpha diversity calculation or relative abundance analyses (McMurdie & Holmes, 2014; Willis, 2019). In case alpha diversity estimates were affected by library size, the alpha diversity analyses were rerun with a data set rarefied to 10,220 reads (the minimum library size after removing small libraries) and there was no change in the results.

| Data sets
The analyses detailed below were conducted on one of two subsets of the data, depending on the questions being explored. The first subset, referred to as the all birds subset, contained birds of all developmental ages (8 day-old nestlings [D8]; 15 day-old nestlings [D15]; adults), and was used to examine the relationship between the gut microbiota, developmental age, and the other life history and environmental variables. The second subset of data contained adults only, and focussed on sex and its interaction with other life history and environmental variables. All analyses were conducted in R version 3.6.3 (R Core Team, 2020).

| Alpha diversity
Linear mixed models (LMM) were used to test the effect of host and environmental factors on alpha diversity. Shannon diversity (Shannon, 1948) measured richness weighted by abundance (the evenness of a community), and Chao1 (Chao, 1984) measured richness, specifically estimating taxa abundance and rare taxa missed from under sampling. Models were fit using the lme4 package (Bates et al., 2015) on each data subset. Significance was determined using Satterwaite's degrees of freedom method (Satterthwaite, 1946) implemented using the lmerTest package (Kuznetsova et al., 2017). Initially the distribution of each alpha diversity metric was assessed graphically and, when obviously skewed, was transformed towards normality as appropriate. Residuals from all models were assessed using DHARMa In the all birds global models, the following variables were included as main effects: age (D8, D15, adult); habitat-type (coniferous, mixed-deciduous); distance to edge, that is, distance from the nest to the nearest woodland edge; maximum brood size of nest; first-egg lay date of nest. All pairwise interactions between these main effects were included, except for the lay date × brood size which was not expected a priori to be related to the microbiota. Sex was excluded because it was unavailable for nestlings.
Mature and immature adult individuals were grouped into a single age bracket "adult" because of convergence issues in the all birds models but not the adults only models. The adults only global models were similar to the all birds models with the addition of a sex variable, and its pairwise interactions with the other variables, and the removal of any habitat or age interactions, due to an imbalance in these factors in the data set (Table S1). Backwards difference coding was used, instead of the default dummy coding, when fitting the age variable which allowed us to compare each age group to the previous age group sequentially, rather than to a single reference level. This contrast scheme gives sum contrasts, so coefficients reflect main effects rather than marginal effects.
Woodland site, nest ID and individual ID were fit to the model as nested random intercepts (in the form woodland site/nest ID/individual ID), to control for nonindependence. Bird ID was dropped from the all birds Chao1 model and woodland site dropped from both adult diversity models due to singular fit warnings, though this had no effect on the model estimates.

| Phylum-level-relative-abundance-GLMMs
The two most abundant phyla (Proteobacteria, mean ± SE: 41.6% ± 2.1%; Firmicutes, mean ± SE: 36.6% ± 2.2%) were modelled separately to determine which variables correlated with changes in the phyla relative abundance, in order to develop a broad sense of how the microbial community changed with varying ecological factors. Binomial generalized linear mixed models (GLMMs) were used, where the binomial response variable was "phyla abundance" per sample. Each sequence read in a sample was scored a "1" if the sequence was from the phylum in question, or scored a "0" if the sequence was from a different phylum. The "phyla abundance" was therefore the summed score divided and weighted by the total number of sequence reads per sample to account for different library sizes. These data were also subset into all birds and adults only subsets, and we used the same model formulas as the alpha diversity models above. The all birds global models did not converge so simplified models were refit only including pairwise interactions with age.
The age variable was coded using backwards difference coding, as in the alpha diversity models.

| Model averaging
Model averaging was performed on the alpha diversity models and the phylum level relative abundance models using the MumIn package (Bartoń, 2020), as outlined by Grueber et al. (2011). Initially full models were fit using lmer or glmer, variables were rescaled using the standardize function in the arm package (Gelman et al., 2021), centring all regression inputs on zero and dividing by 2SDs. Then a model set was generated using the "dredge" command. If multiple models were within the 2AICc cutoff, these were averaged; otherwise the single top model was reported. We also tested model averaging with a less conservative AICc cutoff of 7 (Burnham et al., 2011).
This returned averaged models with more retained terms, but they were all nonsignificant. Moreover, the estimates and significance of the significant terms identified with the 2AICc cutoff did not change, so we only report (conditional) averaged model results from the 2AICc cutoff.

| Model checking and plotting
Model diagnostics for the alpha diversity and phylum level relative abundance models were checked using DHARMa (Hartig, 2019).
DHARMa simulates data based on the model provided and is a better validator than simple residual versus fitted data plots for mixed effects models. Simulated residual plots were made for each of the models in the top model set. In cases where the simulated residuals showed a clear pattern, models were interpreted cautiously. The adults only phyla level models showed overdispersion, which was resolved by refitting the model with an observation level random term following Harrison (2015).
All plots were created using ggplot2 (Wickham, 2016). Alpha diversity and phylum-level results were plotted using the JTools and Interaction (Long, 2019(Long, , 2022 packages, which are based on gg-plot2. The model estimates for the top model in each model set were plotted, not the model averaged results. These plots used partial residuals, which allow the plotting of data accounting for the variables other than the predictor variable of interest and can be better at showing relationships between variables in multivariate mixed models than plotting the raw data, which agreed with the patterns shown in plots of the raw data (not presented).

| Beta-diversity
Beta-diversity measures the dis-similarity between two or more communities, and in our analyses consisted of the pairwise distances between all samples. To address the issue of pseudo-replication in the beta-diversity analyses, due to repeat samples which could not be controlled for using an individual level random effect, the all birds analyses were split into a D8-adult model and a D15-adult model.
For the adults only data set a single model was run. The beta diversity models used the same variables as the equivalent alpha diversity model but only the noninteracting main effects, as these models cannot estimate the main effect and interaction effect at the same time. Beta-diversity was analysed using a compositional approach, which accounts for the proportional nature of high-throughput sequencing data Gloor & Reid, 2016). Low prevalence taxa, that is, those with less than one copy in 5% of samples, were excluded to reduce possible contaminants and sequencing artefacts (Bokulich et al., 2013). The filtered data set was centre-log ratio transformed and then the Euclidean (or Aitchison) distance between samples was calculated (Aitchison, 1983). The beta-diversity of samples was tested using PERMANOVA ("adonis2" function from the vegan package [Oksanen et al., 2019]) after checking variables for homogeneity of dispersion. Models were fit using the "margin" option which tests the marginal effect of each variable while accounting for the other variables in the model. Nest was used as a blocking factor to control for the nonindependence of samples and sequence plate was included as a fixed effect to account for batch effects. Predictor variables were centred and scaled. Some variables in the all birds and adults only models had nonhomogenous dispersions (Table S2 and S3), so significant PERMANOVA results could reflect differences in group variance rather than differences in group means, or could reflect differences in both group variance and group means (Anderson & Walsh, 2013). Beta diversity was visualized using PCA plots of the (pairwise) Euclidian distance between samples, and plotted samples according to their first and second principal component values. Group dispersion was also plotted to visualize differences in group variance, as well as group means.

| Repeatability
Repeatability analyses were conducted on nestlings, for which repeat measures were available across D8 and D15; nestlings with only single measurements were also included to improve estimates (Martin et al., 2011). The rptR package in R (Stoffel et al., 2019) was used to estimate the repeatabilities of individual, nest, and woodland random terms, both on their own, and in combination, to determine the main components of variance underlying microbiota metrics among nestlings. rptR was also used to estimate adjusted repeatabilities, which test whether the components of variance in the random terms were confounded by fixed effects (Nakagawa & Schielzeth, 2010). Thus, repeatabilities for Shannon diversity, Chao1 diversity, phyla-level relative abundance and for each of the first two principal components of beta diversity (using the Euclidian distance matrix as described above), were estimated in the following ways:  Table S4, for phyla level in Table S5, and for beta diversity in Table S6.

| Alpha diversity
There was no relationship between Shannon diversity and age (Table S7a). Chao1 diversity was lower in D15 than in D8 nestlings (−0.35 ± 0.15, p = .02; Table S7b), and there was no difference between D15 nestlings and adult birds (−0.02 ± 0.18, p = .91; Table S7b, Figure 1). The effect of habitat on alpha diversity depended on age -the decrease in Chao1 diversity from D8 to D15 was more pronounced in mixed/deciduous than in coniferous habitats (habitat × age: −0.73 ± 0.32, p = .025; Table S7b, Figure 1). As nest-box distance from the woodland edge increased, so did both Shannon and Chao1 diversity, regardless of age (Shannon: 0.21 ± 0.1, p = .035; Chao1: 0.45 ± 0.2, p = .023; Table S7, Figure 2). Neither brood size nor lay date predicted alpha diversity (Table S7). None of the variables (age, sex, habitat, distance to edge, brood size or lay date) predicted alpha diversity when adults were analysed separately (Table S8).

| Beta diversity
Beta diversity differed across age (Table 1a,b) where community composition became more similar among adult individuals compared with nestlings (Figure 3a,d). In the D8-Adult analysis dispersion tests and plots indicate this age result may be due to a difference in group F I G U R E 1 Alpha diversity varies across life-stages and between habitats. Partial residuals plot for Chao1 diversity regressed on age across habitats (all birds model), means with 95% confidence intervals, shaded circles denote individual samples (n = 246). Red bracket indicates significant difference between D8 and D15 means. Black lines and bracket indicate significant difference between slopes across age and habitat (significant interaction effect). Beta diversity varied with habitat, only when D8 nestlings were included in the analysis (Table 1a-c). Beta diversity varied with lay date and brood size (Table 1a,b), and there was also a tendency for beta diversity to vary with distance to edge (Table 1b), possibly due to a difference in dispersions rather than means (p = < .001, Table S2).
When adults were analysed separately, there was no difference in beta diversity between mature and immature birds or between sexes (Table 1c), and brood size was the only factor to predict beta diversity (Table 1).

| Phylum level relative abundance
D8 nestlings had higher Proteobacteria abundance than D15 nestlings, who in turn had lower Proteobacteria abundance than adults (Table 2a and Figure 4a). The main effect of age also interacted with habitat. In both habitats, Proteobacteria decreased with nestling age (i.e., between D8 and D15), but this effect was stronger for birds developing in mixed/deciduous compared to coniferous habitats, mirroring the pattern above for alpha diversity (Table 2; Figure 4a,b).
Adult birds had greater abundances of Proteobacteria than D15 birds, and this effect was especially pronounced in mixed/deciduous woodlands compared to conifer woodlands (Table 2b; Figure 4a). We note that although the age-habitat interactions are supported by the model, they are not strongly supported by plots (Figure 4a,b) and should be interpreted with caution.
Brood size was negatively related to the abundance of Proteobacteria and positively to the abundance of Firmicutes in nestlings, especially for D15 nestlings, but there was no relationship with either phylum for adults (Table 2a,b; Figure 4c,d). There was no effect of distance to edge on phyla level abundance. Proteobacteria abundance was negatively related to lay date in D8 nestlings, and positively in D15 nestlings and in adults. All of the patterns described for Proteobacteria were found for Firmicutes in the inverse directions ( Figure 4; Table 2). Only lay date predicted phyla level relative abundance in the adults only models. Adult's Proteobacteria increased with lay date (1.5 ± 0.5, p = .006), but did not relate to Firmicutes (−0.7 ± 0.6, p = .223; Table S9).

| Nestling, nest and woodland site repeatability
There were no consistent (repeatable) differences in nestling microbiota among sites (model 1's,  Table S10), indicating that the majority of microbiota variability in nestlings was driven by nest identity, not individual differences. The addition of woodland site as a random effect suggested that up to a quarter of the nest site repeatability may have been driven by differences among sites (Table S10), though the credibility intervals for woodland site repeatability overlapped zero in all six model 6's. Figure 5 includes repeatability estimates and confidence intervals for all six model 6's (i.e., models including site, nest and bird ID, controlling for fixed effects).
Given these repeatability results, we decided to carry out post hoc analyses of how parent traits correlated with offspring traits.
We found that nestling microbiota diversity was predicted only by the mother's Shannon diversity (Table S11a-  Shannon diversity negatively predicted offspring Shannon diversity ( Figure S1). None of the mother's Chao1, Proteobacteria abundance or Firmicutes abundance, or any of the father's microbiota traits predicted the corresponding nestling trait (Table S11e-h). However a post hoc analysis of mean Euclidian distance of each mother from all offspring versus mean distance of each mother from all nonoffspring did not show any difference between community compositions (paired t-test, t = −0.649, p = .52; Figure S3).

| DISCUSS ION
We found that environmental links with the gut microbiota were dependent on age, suggesting differential sensitivity of gut microbiota across life stages in response to habitat type, local environmental gradients and life history traits (Table 3), although experimental perturbations would be necessary to definitively draw this conclusion. The finding that the gut microbiota community stabilizes later in life is consistent with human observational studies (Claesson et al., 2011;Lozupone et al., 2012). Notably, life history, but not environmental effects predicted gut microbiota variation in adults; whereas nestling gut microbiota were predicted by a range of both environmental and life history variations (Table 3). Repeated measures showed that although nestlings appeared to show consistent individual differences in gut microbiota diversity and phylum-level abundance, this was explained entirely by the nest in which they developed, though there was some evidence for site level effects too.
The post hoc parent-offspring analyses found the gut microbiota traits of nestlings were not typically correlated with their parents' traits, although nestlings' Shannon diversity was negatively predicted by their mother's Shannon diversity.

| Environment and age
Only nestling's, not adult's, gut microbiota varied in response to woodland type (Table 3). This is not surprising given microbial community assembly theory (Costello et al., 2012;Coyte et al., 2021) where source microbes vary according to the external environment (e.g., diet, nesting material, habitat) (Chen et al., 2020;Goodenough et al., 2016;Goossens et al., 2022;Koenig et al., 2010) and compete for available niches as the host physiological environment develops (Costello et al., 2012). Given that nestlings immune and digestive systems are probably still developing (Caviedes-Vidal & Karasov, 2001;Stambaugh et al., 2011) and their microbiota is still being established, it is likely that the early microbiota contain transient environmental strains that may be replaced by more stable resident bacteria as the birds mature, similar to the distinct stages found in human infant microbiota development (Stewart et al., 2018). The more variable gut community in nestlings compared to adults is apparent in the alternative beta-diversity plots ( Figure S2) indicating a strong sensitivity of the microbiota to the immediate environment. Due to low recruitment in our study population it was not possible to investigate whether nest effects persisted into maturity.
The alpha diversity of the microbiota was higher for birds in nests further from the edge of the woodland (Figure 2), with a trend for the strength of this association to diminish with age (Table S7b), which we hypothesise could be due to differences in the exposure to microbes via nesting materials (Goodenough et al., 2016) and diet (Davidson, Wiley, et al., 2020) during a time when the microbiota are particularly sensitive to environmental variation. Similarly, D8 nestlings were especially sensitive to the woodland habitat (mixed/ deciduous or coniferous), whereas the effect of woodland on alpha diversity was not present in adults, nor did it differ between adults and D15 nestlings, suggesting that the effects of habitat variation on alpha diversity stabilized early in development. However, for phylumlevel abundance, this stability was not observed until adulthood.
Caterpillars are the main food source for nestling great tits and differ in overall abundance, as well as peak abundance, between deciduous and coniferous woodlands (Balen, 1973;Massa et al., 2004), which may encourage great tits in coniferous woodlands to seek a wider variety of prey (Massa et al., 2004). Future studies should collect microbiota community data and with contemporary DNA metabarcoding data across environmental gradients to determine whether the microbiota covaries with diet across habitats.

| Life history and age
Lay date was positively correlated to Proteobacteria in D15 nestlings and in adults, but negatively in D8 nestlings, while the reverse was true for Firmicutes abundance. Lay date also correlated with variation in community composition (beta diversity) when nestlings and adults were analysed together. The abundance of caterpillars, the primary food source for developing great tits, peaks mid breeding season, the precise timing of which could affect the gut microbiota directly or indirectly through the effects of food availability on stress, immunity or growth (Hooper et al., 2012;Potti et al., 2002;Stothart et al., 2019). Indeed, nestlings with an early or late lay-date may experience stress from lower food availability and higher predation pressure (Naef-Daenzer et al., 2001), and experimentally manipulated glucocorticoids have been reported to affect nestling bird gut microbiota diversity (Noguera et al., 2018). Furthermore, the effect of lay date on stress may depend on nestling age because of the increasing nutritional demands of older nestlings, so that young nestlings from early-nesting parents could face similar conditions to old nestlings from late-nesting parents. If, hypothetically, the positive correlation observed between adult Proteobacteria and lay date was driven by stress, this could explain why the effect of lay date on phylum-level abundance was in opposing directions for D8 and D15 nestlings: D8 nestlings from early nests may have been equally stressed to D15 nestlings from late nests (Figure 4g).
Brood size also correlated with microbiota variation, which we suggest is probably due to larger brood sizes leading to higher stress as a consequence of increased sibling competition (Neuenschwander et al., 2003) and lower food availability (Smith et al., 1988). The effect of brood size on microbiota was especially pronounced for D15 birds, when demand for food is highest and sibling-sibling conflict is greatest. It is possible that larger brood sizes simply have a different reservoir of bacteria due to higher levels of provisioning.
Cross fostering experiments that manipulate the size of broods and compares broods with nestlings from multiple source nests against broods with nestlings from a single source, would test the relative importance of these reservoir and competition hypotheses. Just one other study, to our knowledge, has tested the link between brood size and the gut microbiota in great tits but found no evidence for an effect (Liukkonen et al., 2022). The disparity between these results probably shows how relationships between host and environmental variables and the microbiota may change in different settings or contexts, which will not be apparent in studies taking place during a single season over a relatively small geographic area.
Contrary to our predictions, in the adults only data set, we did not detect any effects of lay date on gut microbiota, nor any sexdependent interactions, despite lay date being strongly determined by female condition (Brown & Brown, 1999). Lay date could also potentially affect the adult gut microbiota indirectly, via increased glucocorticoids related to provisioning effort (Bonier et al., 2011), because provisioning the young is more challenging outside of the peak of food abundance. Although stress hormones have been reported to alter gut microbiota in reproductive female lizards (Sceloporus undulatus), these effects were only observed at specific reproductive stages (MacLeod et al., 2022). In the current study, adults were sampled for gut microbiota near the end of their reproductive cycle (i.e., when chicks were near fledging), and therefore it is possible we missed key reproductive time points where gut microbiota were sensitive to physiological influences associated with lay date. Also, this study uses data from just one breeding season and the relationships TA B L E 2 Phylum level, binomial model results for the all birds (nestlings and adults) subset. Only a single model was retained in the top model sets (<2AICc) so results are not model averaged. The reference level for habitat is "coniferous" which is compared here to "mixed/ deciduous". The age variable is backwards difference coded (see methods for explanation). found here may vary annually depending on environmental conditions. Sampling across multiple years should therefore be an avenue of future research. Nevertheless, brood size, which may also reflect parental quality (Pettifor et al., 2001) and/or higher glucocorticoid demands associated with provisioning (Bonier et al., 2011), predicted gut microbiota community composition in adults of both sexes. Overall, these findings highlight a potentially exciting area of research on reproductive-dependent physiology and parental care, and its interaction with the function of the gut microbiota. Brood size manipulation in conjunction with glucocorticoid measurements would determine whether the effect of brood size on the microbiota was related to parental quality or glucocorticoid levels.

| Repeatability
Our repeatability analyses examined components of gut microbiota variation at individual, local environment and landscape F I G U R E 4 Phyla relative abundance changes with multiple interacting host and environmental factors. Shaded circles denote individual samples (n = 246). Partial residual plots for relative abundance of the two most abundant phyla. Proteobacteria (left column) and Firmicutes (right column) for the all birds phyla-level models showing means with 95% confidence intervals. Point size reflects the weight of that observation, that is, the "number of sequence reads" of the sample. Some points may appear out of range due to "jitter" function which separates overlapping points.  (Mallott, 2022 -Méndez et al., 2022, preprint;Grosser et al., 2019). However we did not see this in our data, lending support to the dominant influence of local environmental factors in nestlings , which we suggest is most probably linked to local food availability. In line with this interpretation, experimental disruption of the microbiota in nestling great tits and blue tits did not have any lasting effect (Diez-Méndez et al., 2022, preprint); the authors provide indirect evidence to suggest that this was most probably due to the parents replenishing the microbiota through feeding and continuous recolonisation from the nest environment. We note that our results do not preclude microbial transfer between parents and nestlings, just that the community and diversity metrics measured here are not comparable between parents and offspring. While we did not detect an effect of host genetics on the nestlings microbiota it is possible that these effects may only become evident at lower taxonomic resolutions or with a larger sample of individual birds. Furthermore these results concern nestlings only and may not hold in adult birds, when the immune system is fully developed-immunogenetic variation has been shown to shape the microbiota of other wild species (Davies et al., 2022;Montero et al., 2021).
Despite lack of evidence for direct transfer of microbiota from the parents themselves (rather than through their provisioned food) to their offspring, we did find a weak negative correlation for the mother's Shannon diversity and that of her offspring. This raises the possibility of a maternal effect on the gut microbiota, which has been demonstrated in wild North American red squirrels (Tamiasciurus hudsonicus) (Ren et al., 2017) and barn swallows (Kreisinger et al., 2017).
One possibility for the negative association is that mothers with high gut microbiota diversity prime their offspring to minimize gut microbiota diversity, akin to parental infections priming offspring to have increased immunity against parasites (Lorenz & Koella, 2011). Our previous study found a negative correlation between gut microbiota diversity and nestling's weight gain, a key proxy for fitness (Davidson et al., 2021), suggesting higher microbiota diversity presents a fitness disadvantage to developing birds. In which case we may expect a selective pressure on mothers to ensure their young do not experience deleteriously high microbial diversity. Maternal investments in egg immunity is also influenced by the bacterial load and diversity of the local environment, affecting nestling development (Jacob et al., 2015;van Veelen et al., 2022).

| CON CLUS IONS
Studies of environmental associations with the gut microbiota in natural systems typically focus on a single ecological feature (e.g., habitat type, season), when in fact multiple interacting ecological factors probably shape the gut microbiota simultaneously. Our results point to links between different aspects of the gut microbiota and multiple aspects of the environment and life history variation: the local nest, sibling interactions, location within woodland plots, the woodland plots themselves, and potentially changes in food availability during the season. We found that environmental effects were important in structuring the gut microbiota, confirming previous work but provide novel evidence for the relative influence of the environment being greater in nestling birds than in adult birds, while life-history traits such as lay date and reproductive effort (brood size) were more important for adults. These effects are to be expected since many of the factors involved are fundamental drivers of ecological processes within populations. Notably, it was the natal as opposed to the adult gut microbiota that was most sensitive to ecological variation, and there was little evidence that the parent microbiota influenced that of their offspring, pointing to a high degree of flexibility during development. Unpicking the mechanisms, causes and consequences of gut microbiota variation among wild populations is an exciting avenue of future research that remains largely unexplored.

ACK N O WLE D G E M ENTS
We thank Iván de la Hera, for assistance with fieldwork, and Niamh Wiley for her help with DNA extractions and sequencing. We also

CO N FLI C T O F I NTER E S T S TATEM ENT
The authors declared no conflict of interest for this article.