Inbred Mouse Populations Exhibit Intergenerational Changes in Intestinal Microbiota Composition and Function Following Introduction to a Facility

Inbred mice are used to investigate many aspects of human physiology, including susceptibility to disease and response to therapies. Despite increasing evidence that the composition and function of the murine intestinal microbiota can substantially influence a broad range of experimental outcomes, relatively little is known about microbiome dynamics within experimental mouse populations. We investigated changes in the intestinal microbiome between C57BL/6J mice spanning six generations (assessed at generations 1, 2, 3, and 6), following their introduction to a stringently controlled facility. Fecal microbiota composition and function were assessed by 16S rRNA gene amplicon sequencing and liquid chromatography mass spectrometry, respectively. Significant divergence of the intestinal microbiota between founder and second generation mice, as well as continuing inter-generational variance, was observed. Bacterial taxa whose relative abundance changed significantly through time included Akkermansia, Turicibacter, and Bifidobacterium (p < 0.05), all of which are recognized as having the potential to substantially influence host physiology. Shifts in microbiota composition were mirrored by corresponding differences in the fecal metabolome (r = 0.57, p = 0.0001), with notable differences in levels of tryptophan pathway metabolites and amino acids, including glutamine, glutamate and aspartate. We related the magnitude of changes in the intestinal microbiota and metabolome characteristics during acclimation to those observed between populations housed in separate facilities, which differed in regards to husbandry, barrier conditions and dietary intake. The microbiome variance reported here has implications for experimental reproducibility, and as a consequence, experimental design and the interpretation of research outcomes across wide range of contexts.


Supplementary Figures and Tables
. Breeding scheme of C57BL/6J inbred mice generations from G1 to G6.   . Relative abundances (proportion) of bacterial taxa that contribute (up to 70% cumulative variation) to the differences observed among the mice generations based on SIMPER analysis. Statistical analyses for each bacterial taxa was performed using ANOVA (litters nested within generation) with Bonferroni correction for multiple comparison between mice generations. Figure S5. Alpha diversity analysis of C57BL/6J inbred mice faecal microbiota from two different facilities, AF1 and AF2, based on (A) microbial richness, (B) Shannon-Wiener index H for microbial diversity and (C) Simpson index (1-D) for microbial evenness, which takes into account microbial richness. (D) The dipersion of microbial community within the two mice populations were also assessed using PERMDISP analysis. Statistical analysis between the groups were performed using the Mann-Whitney test at a significance level of P< 0.05.         Amplicons of the V4 hypervariable region of the bacterial 16S rRNA gene was amplified from faecal DNA extracts as described previously (Choo et al., 2015). To assess LC-MS system variability throughout an analysis run, a pooled biological quality control (PBQC) was prepared by combining 10 µL from each sample.
Samples were analysed in random order interspersed with PBQC measurements every fourth analysis.
Sample (12.5 µg faecal weight equiv) was loaded onto the analytical column in 98% mobile phase A (0.1% aqueous formic acid v/v) and 2% mobile phase B (0.1% formic acid in acetonitrile v/v) at a flow rate of 400 µL/min. The loading conditions were maintained for 0.5 min after which a 12 min linear solvent gradient was applied with a final concentration of 2% A. This concentration was held for 4 min to wash the column after which the solvent mix was returned to starting conditions (98% A) for 2.8 min to re-equilibrate the column ready for the next injection. The same gradient was used for both positive and negative ionisation mode experiments.
Two mass spectrometry experiments were performed, one in negative and one in positive ionisation.
For both modes the instrument was calibrated over the acquisition m/z range of 50 -1500 prior to analysis and mass accuracy was maintained during data acquisition by infusion of a reference solution through the instrument's lockmass channel (200 fmol/µL leucine enkephalin (Sigma-Aldrich, Castle Hill, NSW, Australia) in 1:1 methanol:0.1% aqueous formic acid v/v) which was sampled every 10 s.
Source conditions were optimised for each ionisation mode.

Statistical analysis
Microbial data were analysed for alpha diversity measures (taxa richness, S; Shannon-Wiener index, H; Simpson diversity index, 1-D) of microbial community were determined using PAST (v.3.04) (Hammer et al., 2001). Operational taxonomic unit (OTU) relative abundance was imported into the Primer-E software (v.6, PRIMER-E Ltd, Plymouth, UK) for beta diversity analysis. Bray Curtis similarities were calculated based on the square root-transformed OTU relative abundances, and were used in the non-metric multidimensional scaling (NMDS) ordination plot. Permutational analysis of variance (PERMANOVA) model was used for testing the null hypothesis of no difference (Anderson and Walsh, 2013), based on the parameters permutation of residuals under a reduced model and a type III sum of squares. The bacterial taxa that contributed to the dissimilarities between the groups were determined by SIMPER analysis (Clarke, 1993). PERMDISP was used to assess the dispersion of the microbial community within the groups (Anderson and Walsh, 2013). Comparisons between the microbial and metabolome abundance data was performed using the RELATE analysis (Clarke and Warwick, 2001). Differences in the relative abundance of phyla and genera between generation groups were, where possible, tested for statistical significance based on a nested ANOVA analysis (litters nested in generation), using a type III sum of squares approach on log transformed values in SPSS (v22.0). Multiple pairwise comparisons between groups were performed on the estimated marginal means, with Bonferroni correction applied. Comparisons between generation and within-individual variance was performed using the Kruskal-Wallis test with Dunn's multiple comparison test using GraphPad PRISM 6 (GraphPad Software Inc., La Jolla, USA). Heatmap was generated in R using the ggplots2 (v2.0.0) package (Team, 2015;Wickham and Chang, 2015). The unassigned taxa comprised of less than 0.1% of the relative abundance and were not included in the phyla relative abundance comparison. The DNA and metabolome samples of two mice from G1 and one mouse from G6 failed quality control thresholds and were excluded from microbiota and metabolome analysis, respectively, as the number of observed taxa were very low (less than 20) and the sample was too dilute. G1 mice, and five mice of G6, were not included in assessment of litter effects as data for the maternal origins of these mice were not available.
LC-MS-based metabolomics data were processed using the Progenesis QI software (v2.2, Nonlinear Dynamics, Newcastle Upon Tyne, UK) using default settings. All LC-MS datasets were aligned on retention time, and m/z scales and peak areas were compared. For positive ion data, adduct ions (H + , Na + and K + ) and neutral loss of water ions were identified automatically based on retention time and mass accuracy and combined. For reporting purposes data were filtered to only include markers based on the parameters of ANOVA p-value ≤ 0.005, and at least a 5-fold change in lowest and highest mean intensity between sample groups. Identified markers were exported to EZinfo (v3.0, Umetrics, Umea, Sweden) for OPLS-DA (Orthogonal Projection to Latent Structure-Discriminant Analysis) to generate the S-plots using centred, pareto-scaled data for the identification of discriminant markers between groups. The total variation explained by model (R 2 Y) and predictability of the OPLS-DA model (Q 2 ) were determined according to EZinfo. Significant markers were searched in the METLIN database for potential biomarkers that were within the 10 ppm range of the m/z (or neutral mass if available), and  (Smith et al., 2005). Positive markers with a neutral mass value were only compared against neutral mass searches in METLIN. All statistical analysis were performed at a significance level of 0.05.