1. Natal body mass is a key predictor of viability and fitness in many animals. While variation in body mass and therefore viability of juveniles may be explained by genetic and environmental factors, emerging evidence points to the gut microbiota as an important factor influencing host health. The gut microbiota is known to change during development, but it remains unclear whether the microbiome predicts fitness, and if it does, at which developmental stage it affects fitness traits.
2. We collected data on two traits associated with fitness in wild nestling great tits (Parus major): weight and survival to fledging. We characterised the gut microbiome using 16S rRNA sequencing from nestling faeces and investigated temporal associations between the gut microbiome and fitness traits across development at day 8 (D8) and day 15 (D15) post-hatching. We also explored whether particular microbial taxa were ‘indicator species’ that reflected whether nestlings survived or not.
3. There was no link between mass and microbial diversity on D8 or D15. However, we detected a time-lagged relationship whereby the microbial diversity at D8 was negatively associated with weight at D15, while controlling for weight at D8. Indicator species analysis revealed that while several taxa were unique to birds that either survived or did not survive, there were no universal taxa that were consistently found across all birds within either survival group. This suggests that the presence of particular bacterial taxa may be sufficient, but not necessary for determining future survival, perhaps owing to functional overlap in microbiota.
4. We highlight that measuring microbiome-fitness relationships at just one time point may be misleading, especially early in life. Instead, microbial-host fitness effects should be investigated longitudinally as there may be critical development windows in which key microbiota are established and prime host traits associated with nestling weight. Pinpointing which features of the gut microbial community impact on host fitness, and when during development this occurs, will shed light on population level processes and has the potential to support conservation.

Field monitoring and microbiome sampling

Birds were sampled from nine nest box populations across Co. Cork, Ireland, five of these were mixed/deciduous habitats and four were coniferous habitats (see O'Shea et al. (2018)). We collected 204 faecal samples from 150 nestling great tits from 54 nests (see below for the number of samples that were successfully sequenced and sample size per developmental stage) for 16S rRNA gene sequencing in order to map the gastrointestinal microbial communities. During April-June 2016, nest boxes were monitored to determine lay dates, hatching dates and nestling survival. Nestlings aged 8 days (+/- 1 day) (D8) and 15 days (D15) were placed into sterile holding bags inside a heated holding case. Coffee filters were used to line the bags in order to soak up uric acid from the faeces, as uric acid has the potential to affect downstream sequencing (Khan et al., 1991). Samples were collected from nestlings that defaecated naturally within 15-20 minutes of being out of the nest before being returned to the nest. Although the aim was to have repeated samples for all individuals, not all birds survived, and not all birds produced samples within this time limit. Faecal sacks were opened using a sterile inoculation loop to release the faecal matter and place in a microcentrifuge tube containing 500uL of 100% ethanol. Samples were stored at -20^oC within 8 hours of collection until DNA extraction. D8 birds were weighed, and uniquely individually identified by clipping the tip of one of their nails, avoiding the blood vessel (commonly known as the “quick”). D15 nestlings were weighed and ringed with a unique identifiable metal ring (British Trust for Ornithology) and matched against their unique nail clipping. Nestlings for which we had repeated samples (i.e. for both D8 and D15) are referred to as ‘repeat samples’.

DNA Extractions

Briefly, DNA was extracted from the dried faecal contents of all birds using the Qiagen QIAamp DNA Stool Kit, following the "Isolation of DNA from Stool for Pathogen Detection" protocol (June 2012 edition), with modifications described in (Shutt et al., 2020) to accommodate dried avian faeces. A 0.10 - 0.20 g aliquot of each faecal sample was added to the kit, alongside two negative controls which were carried through to sequencing. The negative controls had low sequence reads and were not included in any subsequent analyses. We note that no decontamination methods were used, though they can be beneficial for samples with low biomass (Eisenhofer et al., 2019). Therefore, we recognise that our sequence data may contain some Amplicon Sequence Variants (ASVs) from the external environment, rather than the gut microbiome, although this should not have systematically biased our results. DNA was stored at -20^oC. Full extraction methods are described in Supporting Information.

Illumina MiSeq sequencing

Full library preparation details are described in Supporting Information and in Davidson et al 2020 Scientific Reports. 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 x 300 cycle kit, following standard Illumina sequencing protocols.

Bioinformatics and statistical analysis

The DADA2 pipeline (Callahan et al., 2016) was used to process the raw sequencing data 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 higher than expected errors. Read errors were estimated before dereplication. Forward and reverse reads were merged to construct 'contig' sequences, these were used to construct a sequence table of ASVs, which 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). This is at 100% sequence identity in contrast to the lower resolution OTU method which groups sequences at 97% identity. ASVs allow greater sensitivity and specificity, better discrimination of ecological patterns than OTU's and are reusable across studies (Callahan et al., 2017).

The DADA2 outputs were assembled into a single Phyloseq object (McMurdie & Holmes, 2013). Sequences identified as mitochondrial or chloroplast were removed. Further filtering of samples and ASVs took place in R. Sample completeness curves were plotted using vegan (Oksanen et al., 2019) and helped determine the lower cut-off for sample reads at 10,000 reads. Low read samples (<10000 reads, 9 samples) were removed leaving 195 (Day-8=81, Day-15=114, repeat samples=41) samples for the analysis. Alpha diversity (both Shannon and Chao1 diversity) was calculated using the ‘estimate_richness’ function from the phyloseq package on the filtered dataset. 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.

See README.txt for description of variables.

Note, some data subsets were centred and scaled before analyses.