Integrated multi-omics analysis identifies microbial and metabolic signatures and drivers of CNS autoimmunity

Abstract

Multiple sclerosis (MS) is an autoimmune disease of the central nervous system (CNS) driven by genetic and environmental determinants. The gut microbiome of people with MS (pwMS) is distinct and influences disease through immunomodulatory metabolite production. Circulating metabolites are altered in pwMS, but identifying microbial-metabolic drivers remains challenging. We previously showed that colonization by the gut commensal Limosilactobacillus reuteri (L. reuteri) exacerbates disease in the experimental autoimmune encephalomyelitis (EAE) model of MS, in a tryptophan-dependent manner. Here, we integrated microbiomic and metabolomic datasets from a longitudinal EAE study utilizing high and low tryptophan diets in mice colonized or not with L. reuteri. Gut microbiome dynamics under short- and long-term alterations in tryptophan bioavailability, were affected by diet, microbiome context, or disease. During short-term dietary intervention, L. reuteri colonization exerted a greater impact on microbiome composition than tryptophan bioavailability. With longer dietary exposure and EAE progression, high dietary tryptophan and L. reuteri colonization synergized to elicit profound microbiota changes, including alterations in Lachnospiraceae, Blautia, and Akkermansia. Integration of metabolomic and microbiomic datasets using joint Robust Aitchison PCA revealed clusters of associated metabolites and microbiota enriched for functional pathways, including bile acid and tryptophan metabolism. Metabolites outperformed microbiota in predicting EAE severity, identifying p-cresols and indoles as top disease-associated metabolites. Treatment with p-cresol or 3-indoleglyoxylic acid exacerbated EAE, enhanced proinflammatory T cell responses, and increased cerebellar pathology. These data demonstrate that dietary responses are shaped by gut microbiome composition and that integrated microbiomic-metabolomic analyses can identify drivers of disease worsening in MS.

Keywords: Limosilactobacillus. reuteri; Multiple sclerosis; cresol; indole; metabolomics; microbiome; multi-omics.

Figure 1.. Schematic of microbiome transplantation, vertical transmission, and dietary tryptophan modulation model, depicting timeline for fecal and serum collection.

Cryopreserved B6 cecal gut microbiome, naturally lacking L. reuteri (B6) or supplemented with 10⁹ CFU L. reuteri (B6 + L. reuteri) was introduced as a one-time inoculation into germ free B6 breeding pairs, which were used to generate experimental offspring that were randomized to high or low tryptophan (Trp) diets. Fecal samples were collected at 1 week post dietary intervention and interrogated by full-length 16S DNA sequencing (FL16S), followed by compositional analyses (see Methods).

Figure 2.. L. reuteri colonization exerts a stronger influence on microbial composition than short-term dietary tryptophan modulation.

Mice colonized by B6 or B6 + L. reuteri microbiomes were randomized to high or low Trp diets, followed by EAE induction (see Fig. 1). Fecal samples were collected at 1 week post dietary intervention (directly prior to EAE induction) and interrogated by full-length 16S DNA sequencing, followed by compositional analyses. (A) Alpha (Shannon), (B) beta (unweighted UniFrac) and (C) unweighted UniFrac PERMANOVA diversity analysis in each gut microbiome × diet treatment group, with significance determined using Wilcoxon rank sum non-parametric (A) or Adonis (B) testing. (D) Stacked bar plots of top 11 most abundant genera as a proportion of total reads. (E) Differentially abundant ASVs (represented as a taxonomic best-hit) between the B6 and B6 + L. reuteri gut microbiomes as determined by LinDA when adjusting for diet, using a cutoff of p_adj ≤ 0.05. Warm colors represent increased abundance in the presence of L. reuteri and cooler colors represent decreased abundance, scaled by the −log(p_adj) × log2 fold-change. Differentially abundant ASVs driven by diet in mice colonized by the B6 (F) or B6 + L. reuteri (G) gut microbiomes, where log2 fold-change indicates increased abundance with a high tryptophan diet when positive and decreased abundance when negative, as determined with LinDA, using a cutoff of p_adj ≤ 0.05. Differentially abundant ASVs in low tryptophan (H) or high tryptophan (I) fed mice driven by microbiome where log fold-change indicates increased abundance within the B6 + L. reuteri gut microbiome when positive and decreased abundance when negative as determined using LinDA with a cutoff of padj ≤ 0.05.

Figure 3.. Baseline microbiota influence the gut microbial compositional response to prolonged changes in tryptophan bioavailability.

Mice colonized by B6 or B6 + L. reuteri microbiomes were randomized to high or low Trp diets, followed by EAE induction (see Fig. 1). Fecal samples were collected at 5 weeks post dietary intervention (30 days post EAE induction) and interrogated by full-length 16S DNA sequencing, followed by compositional analyses. (A) Alpha (Shannon), (B) beta (unweighted UniFrac) and (C) unweighted UniFrac PERMANOVA diversity analyses in each gut microbiome × diet treatment group, with significance determined using Wilcoxon rank sum non-parametric (A) or Adonis testing (B). (D) Stacked bar plots of top 11 most abundant genera as a proportion of total reads. (E) Differentially abundant ASVs (represented as a taxonomic best-hit) between the B6 and B6 + L. reuteri gut microbiomes when adjusted for diet (top), or between low and high tryptophan fed mice when adjusted for gut microbiome (bottom), as determined by LinDA, using a cutoff of p_adj ≤ 0.05. Warm colors represent increased abundance and cooler colors represent decreased abundance, scaled by −log10(padj) × log2 fold-change. Differentially abundant ASVs driven by diet in the B6 (F) or B6 + L. reuteri (G) gut microbiomes, where log fold-change indicates increased abundance with a high tryptophan diet when positive and decreased abundance when negative, as determined using LinDA, with a cutoff of p_adj ≤ 0.05. Differentially abundant ASVs driven by microbiome in low tryptophan (H) or high tryptophan (I) fed mice, where log fold-change indicates increased abundance within the B6 + L. reuteri gut microbiome when positive and decreased abundance when negative as determined using LinDA with a cutoff of padj ≤ 0.05.

Figure 4.. Gut microbiome temporal dynamics are shaped by tryptophan bioavailability and L. reuteri colonization.

Paired microbiome samples collected from mice colonized by B6 or B6 + L. reuteri microbiomes on high and low Trp diets at 1 week and 5 weeks post-dietary intervention were analyzed for the effect of time (see Fig. 1). (A) Alpha (Shannon) and (B) beta (unweighted UniFrac) diversity analysis in each gut microbiome × diet treatment group at the 1-wk and 5-wk post-dietary intervention timepoints. Alpha diversity between timepoints within group is a paired analysis of the absolute change between timepoints using linear models and a t-test for significance. Difference in the magnitude of change in alpha diversity over time between experimental diet-by-microbiome groups represents log-fold change between groups with significance determined by t-test. Beta diversity statistical significance of change over time in each group was evaluated with adonis testing. Differentially abundant ASVs by time in the B6 gut microbiome in mice fed a low (C) or high (D) tryptophan diet, where log fold-change indicates increased abundance at the 5-wk timepoint when positive and decreased abundance when negative, as determined using LinDA with a cutoff of p_adj ≤ 0.05. Differentially abundant ASVs by time in the B6 + L. reuteri gut microbiome in mice fed a low (E) or high (F) tryptophan diet, where log fold-change indicates increased abundance at the 5-wk timepoint when positive and decreased abundance when negative, as determined using LinDA using a cutoff of p_adj ≤ 0.05.

Figure 5.. Integrated microbiomic and metabolomic analysis reveals functionally distinct clusters of microbes and metabolites.

Datasets containing matched fecal microbiome and blood metabolome samples collected at 1 wk and 5 wks post-dietary intervention were integrated using joint robust principal component analysis (jRPCA), which was used to determine sample ordination in principal coordinate space and feature-feature correlation coefficients (see Methods). (A) Emperor biplot of integrated sample clustering in principal coordinate space based on multiomic variation, with arrows indicating the vectors of the top 25 features (microbe ASVs or metabolites) contributing to principal components. (B) Heatmap of correlation coefficients between ASVs (rows) and metabolites (columns), clustered using Euclidean distance and complete linkage, with clusters identified by cutting the trees at fixed heights, resulting in 7 distinct metabolite clusters and 11 ASV clusters. (C) Metabolites in each cluster in (B) were analyzed and annotated by pathway enrichment analysis, with significant enrichment at cutoff of p_adj ≤ 0.05. Networks of correlated ASVs and metabolites within the bile acid (D) or tryptophan (E) clusters. Data was filtered at a threshold of ∣0.99∣ strength of correlation, nodes represent metabolites (grey) or ASVs (with color assignment based on taxonomic class), edge widths represent the strength of correlation, and colors indicate the direction of correlation. (F) ASVs represented in the bile acid or tryptophan networks within each taxonomic class as a percentage of all ASVs in a given class in the total dataset.

Figure 6.. Bile acids, tryptophan, and associated microbiota delineate effects of diet and time in a microbiome-specific manner.

Microbiome and metabolomic datasets were subset by gut microbiome configuration (B6 or B6 + L. reuteri) and integrated using jRPCA. Emperor biplot of integrated sample clustering based on multiomic variation, in the B6 (A) or B6 + L. reuteri (B) gut microbiomes, where vectors indicate the contribution of top 25? microbial and metabolic features to principal components. Top 50 features that delineate samples along principal component (PC) 1 and 2 in the B6 (C and D) or B6 + L. reuteri (E and F) gut microbiomes, segregating samples by dietary context or time.

Figure 7.. Systemic metabolites are better predictors of disease severity than the composition of the gut microbiome.

(A) X-Y scatter plot reflecting total ASV variables of importance by % increase in mean squared error (MSE) and increase in node purity as determined by a Random Forest regression model trained on total ASV level FL16S abundance data to predict cumulative disease score (CDS), with a scatter plot of actual versus predicated values in (B). Top 30 metabolites as variables of importance by % increase in MSE (C) and increase in node purity (D) determined by a Random Forest regression model trained on total metabolomics data to predict CDS, with a scatter plot of actual versus predicated values shown in (E). Top 30 metabolites as variables of importance ranked by % increase in MSE (F) or increase in node purity (G), determined by a Random Forest regression model trained on joint FL16S and metabolomics data, with a scatter plot of actual versus predicated values shown in (H). Model performance was assessed by correlation between actual and predicted values, root mean square error, and the % explained variance. All models were trained using 1000 trees and optimized for feature selection per split and minimum samples per node.

Figure 8.. Indole and cresol metabolites are sufficient to exacerbate CNS disease pathogenesis, promoting proinflammatory T cell responses.

(A) Schematic of experimental model depicting timeline for metabolite treatment with 5 mM p-cresol, 1 mM 3-indoleglyoxylic acid, (IGoxA), or vehicle control in the drinking water, relative to EAE induction in B6 Jax mice at 8 weeks of age (n=5M/group). EAE disease courses for each metabolite treatment are represented as mean daily clinical score with overall significance determined by Friedman’s parametric two-way ANOVA for (B) classical EAE scoring and (C) AR-EAE scoring, assessing the time*treatment interaction effect. CDS representing the total sum of all daily EAE scores over a 30 day disease course for metabolite treatment groups and is shown in (D), with significance evaluated using Mann-Whitney non-parametric test. (E) % starting weight represents mean daily percentage within each group with overall significance determined by Friedman’s parametric two-way ANOVA, assessing the time*treatment interaction effect. At day 30 post EAE induction, spinal cord-infiltrating leukocytes were isolated and analyzed by flow cytometry. Frequencies of the major leukocyte populations within the CD45⁺ population and total cell counts of (F and G) CD11b⁺, (H and I) CD19⁺, (J and K) TCRβ⁺, (L and M) TCRγδ⁺, (N and O) CD4⁺ T cells, and (P and Q) CD8⁺ T cells. Cytokine production as a frequency of TCRγδ⁺ T cells is shown for IFNγ (R) and IL-17 (S). Significance was determined by one-way ANOVA at p_adj ≤ 0.05.

Figure 9.. Indole and cresol metabolites promote cerebellar brain pathology.

EAE and metabolite treatments were conducted as outlined in Figure 8A. At day 30 post EAE induction, brains were collected and processed for staining with H&E with or without LFB. Histopathologic evaluation of control, p-cresol, and IGoxA (n=4–5M/group) was performed as described in Methods. (A and B) Cerebellum inflammation scores by strain and corresponding H&E representative images. (C and D) Cerebellum demyelination scores by strain and corresponding LFB with H&E representative images. Cerebellum images (B and D) were captured at 5× objective. Scale bar: 200 μm. For all images, the arrows mark regions of inflammatory infiltrates or demyelination. Significance of differences was determined by ordinary 1-way ANOVA, with Fishers LSD multiple-comparison test at p ≤ 0.05.