DNA extraction, 16S rRNA gene sequencing, and gut microbiota analyses

DNA from fecal samples were extracted using a QIAamp Fast DNA Stool Mini Kit (Qiagen Cat# 51604) following the manufacturer’s instructions. Amplicon sequencing of the 16S rDNA V3-V4 region of bacterial communities was performed as previously described,⁵⁶ utilizing 2x250bp paired-end Illumina MiSeq chemistry (Genewiz; Leipzig, Germany). Raw sequencing reads were processed using an in-house 16S rRNA processing pipeline, employing USEARCH (64 bit; v 8.1). Briefly, paired-end reads were merged and filtered using a < 0.5 expected error rate per nucleotide and total length. Reads were dereplicated, and singletons removed, following the trimming of the forward and reverse primers (“-stripleft 17” and “-stripright 21”, respectively). Operational taxonomic units (OTUs) were clustered at 97% identity, and reference-based chimera removal was performed using UCHIME. OTUs were assigned taxonomic information by aligning reads to the RDP Gold database using the RDP Classifier (v 2.12).

The 16S rRNA reads assigned taxonomic information was converted into a count matrix and imported into R Studio (v 3.6.1). The 16S data was analyzed using the metadata accompanying the relevant mouse trials and the RT-qPCR. The 16S rRNA OTU reads aligned per OTU were converted into relative abundances using the “funrar” package.⁵⁷ Dataframes and matrices were manipulated as necessary using the “reshape2” package.⁵⁸ Publication quality images were generated using the “ggplot2” and “ggpubr” packages.^59,60 Boxplots represent the standard Tukey representation, with boxes representing the 25th, 50th (median) and 75th interquartile range (IQR) percentiles and the whiskers encompassing values within 1.5× the IQR. The specific values of each boxplot are overplotted as opaque gray circular points using the ggplot “geom_jitter” function, whereas boxplot outliers are represented as solid square points.

The alpha- and beta-diversity of the samples were calculated using the R packages “vegan” and “phyloseq”.^61,62 All alpha-diversity values presented use Simpson’s index, while beta-diversity separation was performed using Canberra distances with PCoA ordination. A color palette for the composition bar plots was obtained through the “pals” R package.⁶³ The log₂ fold change of epithelial and inflammatory biomarker genes, and 16S rRNA taxa, was performed using the “DESeq2” package after scaling values to integers.⁶⁴ Correlations between microbial taxa and epithelial and inflammatory biomarker genes were performed using base R’s “stats” package,⁶⁵ with plots generated using the “corrplot” package.⁶⁶ Individual images were manipulated into their final multipanel display using Inkscape (v 1.1.2).