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.