Beta diversity analysis

The unrarefied reads (filtered to remove samples with < 10,000 reads) were used. Samples were processed to remove exceptionally rare taxa that could disproportionately influence beta diversity metrics. As such, ASVs were excluded if they had fewer than 50 reads in total across all samples and/or were present in less than 2% of samples. This retained 3,057 ASVs across 586 samples. ASV abundances were then transformed using the Centered Log Ratio (CLR) transform function in microbiome 1.12.0 [75]. The CLR transformation produces values that are scale invariant (i.e. not influenced by differences in library sizes across samples) and controls for the compositional nature of microbiome datasets [76]. The consistency of beta diversity estimates across different extractions and sequencing runs was assessed as described for Shannon diversity, but using pairwise Euclidean distances calculated using the CLR transformed ASV abundances.

To identify whether compositional differences in the GM were associated with host age class, sex, or environmental variables, sequencing and catch duplicates were first removed (as described above). Samples taken from nestlings and floaters were also removed (as described above), leaving 450 samples from 309 individuals. A Euclidean distance matrix was then calculated using the CLR-transformed ASV abundances. To quantify differences in beta diversity between groups of samples, a Permutational Analysis of Variance (PERMANOVA) was performed using the adonis2 function within the R package vegan 2.5.7 [77, 78], with 9999 permutations. Individual age class (fledgling, old fledgling, sub-adult, adult) and sex were included as predictors, as well as territory quality (as a proxy for dietary differences across the island and local competition for food) and the corresponding sampling field period (Major 2017, Minor 2018, Major 2018, Minor 2019, Major 2019, Minor 2020). Bird ID was included as a blocking factor to control for repeated sampling. The function betadisper was used to check for the homogeneity of group dispersion values [77, 78]. Pairwise PERMANOVA analyses were conducted using pairwiseAdonis 0.0.1 [79]—the Benjamini and Hochberg method [66] was used to correct P values for multiple testing as part of this method. Differences in GM composition across predictor variables were visualised via a Principal Components Analysis (PCA). A second analysis was also performed using ASV abundances that were instead transformed using the Phylogenetic Isometric Log Ratio transformation (PhILR) in the R package philr 1.14.0 [80]. This transformation controls for the compositionality of the data but also preserves information about the relatedness of ASVs, thus providing a phylogenetically aware measure of beta diversity [80].

Samples were filtered as above (see alpha diversity analysis), leaving a total of 425 samples from 296 individuals. As GM composition differed across age classes and adult birds were found to be significantly heavier than juvenile individuals (fledgling, old fledgling and sub-adult age classes, Additional file 1: Table S2), subsequent analyses were carried out separately for juveniles (205 samples from 175 individuals) and adults (220 samples from 165 individuals); this was to ensure that any confounding influence of age class and body condition on GM composition could be separated. For juvenile age classes, a PERMANOVA analysis was conducted using a Euclidean distance matrix of CLR-transformed ASV abundances, with 9,999 permutations. Body condition was included as a predictor variable, by extracting the residuals of a regression of body mass on tarsus length, controlling for the time of day, separately for males and females. Individual sex, territory quality and the corresponding sampling field period (Major 2017, Minor 2018, Major 2018, Minor 2019, Major 2019, Minor 2020) were included as additional predictors. Age class (fledgling, old fledgling or sub-adult) was also included in the juvenile model. Bird ID was included as a blocking factor to control for repeated sampling. A second analysis was also performed using ASV abundances that were instead transformed using the PhILR transformation [80].

To identify compositional differences in the GM between individuals that survived versus those that died, sequencing and catch duplicates were removed, as were samples taken from nestlings, floaters and those collected in 2020. The remaining samples were also filtered to retain the latest sample per individual (as described above). Post-filtering, 264 samples containing 2,900 ASVs were retained (226 individuals that survived, 38 individuals that died). As for the analysis of beta diversity and body condition, analyses were carried out separately for juveniles (116 samples/individuals; 17 died, 99 survived) and adults (148 samples/individuals; 21 died, 127 survived), since survival probability in the Seychelles warbler is significantly lower in the first year of life (encompassing fledglings to sub-adults) [44]; as such, separating juveniles and adults ensures that changes in GM structure that are associated with age and survival are not confounded. For each age group (juvenile or adult) a PERMANOVA analysis was performed, using a Euclidean distance matrix of CLR-transformed ASV abundances, with 9,999 permutations. Survival to the next breeding season (yes, no), sex, territory quality and the corresponding sampling field period (Major 2017, Minor 2018, Major 2018, Minor 2019, Major 2019) were included as predictors. Age class (fledgling, old fledgling or sub-adult) was also included in the juvenile model. Differences in GM composition across samples were visualised via PCA. As above, a second analysis was also performed using a Euclidean distance matrix calculated from PhILR-transformed ASV abundances to assess whether communities were phylogenetically distinct across groups.

To establish whether specific ASVs were differentially abundant across groups of individuals, an Analysis of Compositions of Microbiomes with Bias Correction (ANCOM-BC) was carried out, using the ANCOMBC 1.1.5 package in R [81]. Differential abundance was tested between groups of adult individuals that survived, versus those that died, whilst controlling for sex and sampling season. As part of ANCOM-BC, the Benjamini and Hochberg method was used to correct P values for multiple testing [66]. A cut-off of P_adj < 0.05 was used to assess significance.