Gut microbiome diversity, composition, and relative abundance analyses
To account for uneven sequencing depth across samples, samples were rarefied to the depth of our lowest sample (5,691 for the diet manipulation experiment; 1,077 for the diet selection experiment) for alpha and beta diversity analyses. We used the vegan package (Oksanen et al., 2022) to quantify two measures of alpha diversity, the observed number of ASVs (richness) and the Shannon diversity index (which accounts for both richness and evenness). To determine whether the observed microbial richness varied in response to treatment, we ran a generalized linear mixed models (GLMM) fitted with a negative binomial (nbinom1) distribution using the package glmmTMB (Brooks et al. 2017) with treatment, day, and the interaction between treatment*day as fixed effects. For the Shannon diversity index, we ran a GLMM fitted with a gaussian structure and included treatment, day, and the interaction between treatment*day as fixed effects. Because birds were sampled over time, we included bird ID as a random effect in all models. We used the DHARMa package (Hartig 2019) to plot residuals and confirm the suitability of each model and the Anova function in thecar package (Fox and Weisberg 2018) to determine significance. To assess beta diversity, we first subjected the unrarefied dataset to Cumulative Sum Scaling (CSS) normalization (Paulson et al. 2013) using the metagenomeSeq package following the methods outlined in Maraci et al. 2021. To account for compositional variations in each dataset, data were log (x+0.0001)-transformed and later corrected by subtracting the log of the pseudo count (Thorsen et al. 2016). For beta diversity analyses, we computed Bray-Curtis (Bray and Curtis 1957), unweighted UniFrac (Louzupone and Knuight 2007) and weighted UniFrac (Lozupone et al. 2007) distances. To visualize dissimilarities between treatments, we used a principal coordinate analysis (PCoA) using the ordinate function in phyloseq (McMurdie & Holmes, 2013). To assess similarity between treatments and over time for each beta diversity metric, we conducted PERMANOVA tests using the adonis2 function in the veganpackage (Oksanen et al., 2022). Treatment, day sampled, and the interaction between treatment and day were included as fixed effects in all PERMANOVA tests. To account for repeated sampling over time, bird ID was included under the adonis2 strata option. All P-values were adjusted for false discovery rate with a Benjamini–Hochberg correction where significance was determined as Padj < 0.05.
Finally, to investigate whether bacterial taxa differ in abundance across treatments, we calculated the relative abundances of phyla and genera from the unrarefied datasets. For analyses, phyla and genera with mean abundances <1% were lumped into an “Other” category. This data stringency limited analyses to the top four most abundant phyla (Campylobaterota, Firmicutes, Proteobacteria, Actinobacteriota) and the top seven most abundant genera (Campylobacter ,Weissella , Helicobacter , Acinetobacter ,Ligilactobacillus , Pseudomonas , Achromobacter ) in experiment 1 (diet manipulation). In experiment 2 (diet selection), this data stringency limited analyses to the five most abundant phyla (Campylobaterota, Firmicutes, Proteobacteria, Actinobacteriota, Bacteroidota) and the top ten most abundant genera (Campylobacter , Helicobacter , Corynebacterium ,Ureaplasma , Ligilactobacillus , Atopobium ,Catellicoccus , Gallibacterium , Veillonella,Weissella ) in injected birds (LPS, saline), and the top thirteen most abundant genera (Campylobacter , Helicobacter ,Pseudomonas , Serratia , Veillonella,Gallibacterium , Ligilactobacillus , Atopobium ,Corynebacterium , Enterobacter , Sphingobacterium ,Staphylococcus , Weissella ) in focal birds (LPS-focal, saline-focal). To compare abundances of bacterial taxa between treatments and over time, we ran non-parametric Kruskal-Wallis tests using the kruskal.test function in R. P-values were adjusted for false discovery rate with a Benjamini–Hochberg correction where significance was determined as Padj < 0.05.