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

The gut microbiome has important consequences for fitness, yet the complex, interactive nature of ecological factors that influence the gut microbiome has scarcely been investigated in natural populations. We sampled the gut microbiota of wild great tits (Parus major) at different life stages and across multiple conifer and mixed woodland fragments, allowing us to evaluate multiple factors that relate to within-individual gut microbiota acquisition, including habitat type, nest position and life history traits. The gut microbiota varied with both environment and life-history in ways that were largely dependent on age. Notably, it was the nestling, as opposed to the adult gut microbiota that was most sensitive to ecological variation, pointing to a high degree of developmental plasticity. Individual nestling differences in gut microbiota were consistently different (repeatable) from one to two weeks of life, driven entirely by the effect of sharing the same nest. Our findings point to important early developmental windows in which the gut microbiota are most sensitive to environmental variation and suggest reproductive timing, and hence parental quality or food availability, interact with the microbiome.


7
estimating individual repeatability using individual as a random effect in a mixed model 150 framework (Nakagawa & Schielzeth, 2010). Mixed models also enable an examination of 151 whether individual differences are independent of other individual or environmental factors 152 that can also be included as random or fixed effects in the same model. Despite its 153 widespread use in animal personality and cognition literature (Bell et al., 2009;Cauchoix et 154 al., 2018), to our knowledge this approach has yet to be applied to identify the sources of 155 variation in the gut microbiota. The DADA2 pipeline was used to process the raw sequencing data (Callahan et al., 2016) in 227 R version 3.5 (R Core Team, 2019). Sequence quality was visually inspected. Sequences 228 were trimmed to remove adapters and lower quality reads (median quality scores below  30 threshold) at the extremities of the sequence and filtered to remove sequences with 230 expected errors above 1. Read errors were estimated before dereplication. Forward and 231 reverse reads were merged to construct 'contig' sequences, which were used to construct a 232 sequence table of Amplicon Sequence Variants (ASV's), which in turn counts the number of 233 times each unique sequence is detected. The previous steps were performed for each run 234 separately. Then the separate sequence tables were merged and chimeras removed using the 235 'consensus' method. Taxonomy was assigned to each ASV by RDP's Naive Bayes Classifier 236 (Wang et al., 2007) against the Silva reference database (version 132) (Quast et al., 2012). 237 This method groups sequences with 100% sequence identity in contrast to the lower 238 resolution OTU method which groups sequences at 97% identity. ASV's allow greater 239 sensitivity and specificity, better discrimination of ecological patterns than OTU's and are 240 reusable across studies (Callahan et al., 2017). Linear mixed models (LMM) were used to test the effect of host and environmental factors 275 on alpha diversity. Shannon diversity (Shannon, 1948) measured richness weighted by 276 abundance (the evenness of a community), and Chao1 (Chao, 1984) measured richness, 277 specifically estimating taxa abundance and rare taxa missed from under sampling. Models 278 were fit using the lme4 package (Bates et al., 2015) on each data subset. Significance was 279 determined using Satterwaite's degrees of freedom method (Satterthwaite, 1946) 280 implemented using the lmerTest package (Kuznetsova et al., 2017). The distribution of each 281 alpha diversity metric was assessed graphically and transformed towards normality as 282 appropriate. 283

284
The following variables were included as main effects: age; sex; habitat-type: 285 coniferous/mixed; distance to edge, i.e. distance from the nest to the nearest woodland edge; 286 maximum brood size of nest; first-egg lay date of nest. In the all birds global models, all 287 pairwise interactions between these main effects (except for lay date × brood size, which 288 were strongly colinear; sex also excluded because it was unavailable for nestlings) were 289 included. In the adults only global models, all pairwise interactions were included, except for 290 lay date × brood size and any habitat or age interactions, due to an imbalance in these factors 291 in the dataset, (see table S1). Backwards difference coding was used, instead of the default 292 dummy coding, when fitting the age variable which allowed us to compare each age group to 293 the previous age group sequentially, rather than to a single reference level. This contrast 294 scheme gives sum contrasts, so coefficients reflect main effects rather than marginal effects. 295 Woodland site, nest ID and individual ID were fit to model as nested random intercepts (in 296 the form woodland site/nest ID/individual ID), to control for non-independence. Bird ID was 297 dropped from the all birds Chao1 model and woodland site dropped from both adult diversity 298 models due to singular fit warnings.

Phylum-level-Relative-Abundance-GLMMs 301
The two most abundant phyla (Proteobacteria, mean ± SE: 41.6% ± 2.1%; Firmicutes, mean 302 ± SE: 36.6% ± 2.2%) were modelled separately to determine which variables correlated with 303 changes in the phyla relative abundance, in order to develop a broad sense of how the 304 1 5 was calculated (Aitchison, 1983). The beta-diversity of samples was tested using 350 PERMANOVA (adonis2 function from the vegan package (Oksanen et al., 2019)) after 351 checking variables for homogeneity of dispersion. Models were fit using the 'margin' option 352 which tests the marginal effect of each variable while accounting for the other variables in the 353 model. Nest was used as a blocking factor to control for the non-independence of samples 354 and sequence plate was included as a fixed effect to account for batch effects. Predictor 355 variables were centred and scaled. For the all birds dataset two models were run: one with 356 only non-interacting fixed effects, and another model which included the main interaction of 357 interest (age × habitat), as PERMANOVA cannot calculate the marginal effect of fixed 358 effects from the interactions in the model. For the adults only dataset a single model was run 359 with only non-interacting fixed effects. No variable, except habitat, had homogenous 360 dispersion when all birds were considered (table S2), so significant PERMANOVA results 361 could reflect differences in group variance rather than differences in group means, or could 362 reflect differences in both group variance and group means (Anderson & Walsh, 2013). 363 When adult birds were considered separately age, sex and habitat all had homogenous 364 dispersions (table S3). Beta diversity was visualised using PCA plots of the (pairwise) 365 Euclidian distance between samples, and plotted samples according to their first and second 366 principal component values. 367

368
Repeatability 369 The rptR package (Stoffel et al., 2019) was used to decompose the components of variance 370 underlying alpha diversity in nestlings measured across D8 and D15, and to examine the 371 extent to which consistent differences among nestlings were confounded by the fixed effects 372 of age and habitat, and the random effects of nest or woodland site. Specifically, repeatability 373 of Shannon and Chao1 diversity was calculated at the individual and group level in the suggests that the distinction between the ages was most pronounced in the mixed/deciduous 409 woodland (figure S1), though dispersion tests indicate this result could be due to a difference 410 in group dispersions, group means or both (table S2). Beta diversity varied with brood size 411 and distance to edge, and there was also a tendency for beta diversity to vary with lay date 412 (table 1a). When adults were analysed separately, there was no difference in beta diversity 413 between mature and immature birds or between sexes (table 1c), and brood size was the only 414 factor to predict beta diversity (table 1). 415 416 Phylum level relative abundance 417 D8 nestlings had higher Proteobacteria abundance than D15 nestlings, who in turn had lower 418 Proteobacteria abundance than adults (table 2a and figure 4A). The main effect of age also 419 interacted with habitat. In both habitats, Proteobacteria decreased with nestling age (i.e. 420 between D8 and D15), but this effect was stronger for birds developing in mixed/deciduous 421 compared to coniferous habitats, mirroring the pattern above for alpha diversity (table 2,  422 figure 4A, 4B). Adult birds had greater abundances of Proteobacteria than D15 birds, and this 423 effect was especially pronounced in mixed/deciduous woodlands compared to conifer was not observed until adulthood. Future studies should collect microbiota community data 500 and with contemporary DNA metabarcoding data across environmental gradients to 501 determine whether the microbiota covaries with diet across habitats. 502 503 Life history and age 504 Lay date was positively correlated to Proteobacteria in D15 nestlings and in adults, but 505 negatively in D8 nestlings, while the exact reverse was true for Firmicute abundance. The 506 abundance of caterpillars, the primary food source for developing great tits, peaks mid 507 breeding season, the precise timing of which could affect the gut microbiota directly, or 508 indirectly through the effects of food availability on stress, immunity or growth (Hooper et 509 al., 2012;Potti et al., 2002;Stothart et al., 2019). Indeed, nestlings with an early or late lay-510 date may experience stress from lower food availability and higher predation pressure (Naef-511 Daenzer et al., 2001), and experimentally manipulated glucocorticoids have been reported to 512 affect nestling bird gut microbial diversity (Noguera et al., 2018). Furthermore, the effect of 513 lay date on stress may depend on nestling age because of the increasing nutritional demands 514 of older nestlings, so that young nestlings from early-nesting parents could face similar 515 conditions to old nestlings from late-nesting parents. If, hypothetically, the positive 516 correlation observed between adult Proteobacteria and lay date was driven by stress, this 517 could explain why the effect of lay date on phylum-level abundance was in opposing 518 directions for D8 and D15 nestlings: D15 nestlings from early nests may have been equally 519 stressed to D8 nestlings from late nests (figure 4G). 520 521 Brood size also correlated with microbiota variation, which we suggest is likely due to larger 522 brood sizes leading to higher stress as a consequence of increased sibling competition 523 brood size on microbiota was especially pronounced for D15 birds, when demand for food is 525 highest and sibling-sibling conflict greatest. It is possible that larger brood sizes simply have 526 a different reservoir of bacteria due to higher levels of provisioning. Cross fostering 527 experiments that manipulate the size of broods, and compares broods with nestlings from 528 multiple source nests against broods with nestlings from a single source, would test the 529 relative importance of this reservoir and competition hypotheses. 530 531 Contrary to our predictions, in the adults only dataset, we did not detect any effects of lay 532 date on gut microbiota, nor any sex-dependent interactions, despite lay date being strongly 533 determined by female condition (Brown & Brown, 1999). Lay date could also potentially 534 affect the adult gut microbiota indirectly, via increased glucocorticoids related to 535 provisioning effort (Bonier et al., 2011), because provisioning the young is more challenging 536 outside of the peak of food abundance. Although stress hormones have been reported to alter 537 gut microbiota in reproductive female lizards (Sceloporus undulatus), these effects were only 538 observed at specific reproductive stages (MacLeod et al., 2022). In the current study, adults 539 were sampled for gut microbiota near the end of their reproductive cycle (i.e. when chicks 540 were near fledging), and therefore it is possible we missed key reproductive time points 541 where gut microbiota were sensitive to physiological influences associated with lay date. 542 Nevertheless, brood size, which may also reflect parental quality (Pettifor et al., 2001) and/or 543 higher glucocorticoid demands associated with provisioning (Bonier et al., 2011), predicted 544 gut microbiota community composition in adults of both sexes. Overall, these findings 545 highlight a potentially exciting area of research on reproductive-dependent physiology and 546 parental care, and its interaction with the function of the gut microbiome. Brood size 547 manipulation in conjunction with glucocorticoid measurements would determine whether the 548 effect of brood size on the microbiota was related to parental quality or glucocorticoid levels. The repeatability analysis revealed the importance of the immediate nest environment in 552 structuring the microbiota, whereas the larger-scale effect of woodland site was negligible. 553 Human parents provide their offspring with beneficial microbes acquired over their own 554 lifespan and hence provide their offspring a fitness advantage (Sprockett et al., 2018). These 555 microbes prime them for their environment during early life and are particularly influential 556 due to priority effects, where the timing of colonisation influences future species interactions 557 (Sprockett et al., 2018). Notably, individual-level repeatability of microbial diversity was 558 zero when nest identity was controlled, where individual's gut microbiota were not following 559 individual trajectories but were converging on their nestmates. The importance of the nest in 560 determining microbiota diversity and structure has been highlighted in this species and others 561 previously (Campos-Cerda & Bohannan, 2020;Teyssier, Lens, et al., 2018). If genetics were 562 an important driver of the microbiota in nestlings we would expect to see some level of 563 repeatability, even when controlling for nest as siblings are expected to share only 50% of 564 their DNA on average. These results point to considerable plasticity in microbiota assembly 565 during early development, and given the negligible evidence for positive correlation between 566 any parent and offspring microbiome trait, is further support to the greater influence of local 567 environmental factors, most likely linked to local food availability, over intrinsic or host 568 factors, in avian microbiota (Song et al., 2020). Despite these findings, the negative 569 correlation for the mother's Shannon diversity raises the possibility of a maternal effect on 570 the gut microbiota, which has been demonstrated in wild North American red squirrels 571 (Tamiasciurus hudsonicus) using a quantitative genetic approach (Ren et al., 2017). Studies of environmental associations with the gut microbiome in natural systems typically 575 focus on a single ecological feature (e.g. habitat type, season), when in fact multiple 576 interacting ecological factors likely shape the gut microbiome simultaneously. Our results 577 point to links between different aspects of the gut microbiota and multiple aspects of the 578 environment and life history variation: the local nest, sibling interactions, location within 579 woodland plots, the woodland plots themselves, and potentially changes in food availability 580 during the season. These effects are to be expected since many of the factors involved are 581 fundamental drivers of ecological processes within populations, and yet notably, it was the 582 natal, as opposed to the adult gut microbiota that was most sensitive to ecological variation, 583 pointing to a high degree of developmental plasticity. Unpicking the mechanisms, causes and 584 consequences of gut microbiome variation among wild populations is an exciting avenue of 585 638

651
S h a n a h a n , N a t u r e , 5 0 5 (  7  4  8  4  ) ,

834
S e x d  i  f  f  e  r  e  n  c  e  s  a  n  d  h  o  r  m  o  n  a  l  e  f  f  e  c  t  s  o  n  g  u  t  m  i  c  r  o  b  i  o  t  a  c  o  m  p  o  s  i  t  i  o  n  i  n  m  i  c