Multi-omic evaluation of metabolic alterations in multiple sclerosis identifies shifts in aromatic amino acid metabolism.

The circulating metabolome provides unique insights into multiple sclerosis (MS) pathophysiology, but existing studies are relatively small or characterized limited metabolites. We test for differences in the metabolome between people with MS (PwMS; n = 637 samples) and healthy controls (HC; n = 317 samples) and assess the association between metabolomic profiles and disability in PwMS. We then assess whether metabolic differences correlate with changes in cellular gene expression using publicly available scRNA-seq data and whether identified metabolites affect human immune cell function. In PwMS, we identify striking abnormalities in aromatic amino acid (AAA) metabolites (p = 2.77E-18) that are also strongly associated with disability (p = 1.01E-4). Analysis of scRNA-seq data demonstrates altered AAA metabolism in CSF and blood-derived monocyte cell populations in PwMS. Treatment with AAA-derived metabolites in vitro alters monocytic endocytosis and pro-inflammatory cytokine production. We identify shifts in AAA metabolism resulting in the reduced production of immunomodulatory metabolites and increased production of metabotoxins in PwMS.

Introduction

Multiple sclerosis (MS) is an inflammatory and neurodegenerative disorder of the central nervous system (CNS) The etiology of MS is multifactorial and involves multiple levels of biological interactions with genetic and environmental contributors., However, a detailed understanding of these interactions or underlying mechanisms remains rudimentary. Molecular profiling of circulating small molecules using metabolomics integrates many of these systems, as an individual’s metabolic phenotype reflects an intersection between environmental sources of variation such as lifestyle characteristics (e.g., diet), upstream genetic influences, and activity of the gut microbiota., Thus, the assembly of a detailed global map of functional relationships in MS using circulating biologic intermediates may provide a valuable step forward in improving our understanding of potential contributors to MS pathogenesis.

Preliminary studies suggest that metabolic alterations exist in people with MS (PwMS) with respect to global metabolomic profiles, specific pathways, and individual metabolites.6, 7, 8, 9 These results are consistent with observations from metabolomics studies of other neurological disorders, including Alzheimer’s and Parkinson’s diseases, and other autoimmune disorders.10, 11, 12 However, with respect to MS, these studies have been restricted to small patient populations and limited arrays of metabolites. Critically, most prior studies have not evaluated the links between metabolomic differences and patient characteristics, including measures of disease severity or other candidate biological drivers of disease.

In the present study, we compare nearly 1,000 detailed metabolomic profiles from PwMS and healthy people from across the age spectrum using individual metabolite and pathway-level analyses. Our findings highlight distinct abnormalities in aromatic amino acid metabolism that imply an altered balance of immunomodulatory metabolites (e.g., arylhydrocarbon receptor [AhR] and hydroxycarboxylic acid-3 [HCA3] receptor agonists) and a shift toward production of known metabotoxins (e.g., indole acetate, phenylacetylglutamine, p-cresol sulfate/glucuronide). Metabolomic alterations, particularly in these immunomodulatory amino acids and metabotoxins, were associated with disease severity. We then integrated publicly available single-cell transcriptomics data to identify alterations in genes involved in aromatic amino acid metabolism in MS peripheral blood and cerebrospinal fluid (CSF)-derived immune cells. Lastly, we tested the effects of the identified AAA-derived metabolites on human peripheral immune cells and noted changes in the production of pro-inflammatory cytokines, as well as endocytosis in monocytes (as a candidate cell type), suggesting functional roles for metabolites that are altered in MS and are related to disease severity.

Results

After applying several quality-control steps (Figures S1 and S2; STAR Methods), we included 954 metabolomic profiles from 756 individuals with MS or healthy controls (HC) (514 MS patients and 241 HC), in which 269 metabolites were reliably measured. Participants were aged 40.26 ± 15.46 years and were predominantly female (73.49%) and non-Hispanic whites (84.21%; Table 1; Figure S1; Data S1A). With respect to MS, the majority of participants had relapsing-remitting MS (RRMS; 72.71%), disease duration was 13.45 ± 10.63 years, and they were, on average, moderately disabled (median Expanded Disability Status Scale [EDSS] score = 3.0; 27% reported using a cane).

Table 1

Characteristics of included study participants

Disease status	HC	MS
No. samples	317	637
No. of participants	241	515
Age, y, mean (SD)	35.85 (15.71)	42.54 (14.91)
Age range (min, max)	7, 84	7, 76
Female gender, n (%)	224 (70.66)	468 (73.47)
White race, n (%)	253 (79.81)	542 (85.09)
Disease duration, y, mean (SD)	–	13.45 (10.63)
Progressive MS, n (%)	–	205 (32.18)
EDSS (417 scores available), median (IQR)	–	2.50 (1.50, 5.00)
Use of cane, n (%)	–	97 (23.2)
Disease-modifying therapy, n (%)	–	–
No treatment	–	195 (30.61)
Interferon-β	–	60 (9.42)
Glatiramer acetate	–	130 (20.41)
Dimethyl fumarate	–	36 (5.65)
Teriflunomide	–	2 (0.31)
Natalizumab	–	82 (12.87)
Alemtuzumab	–	2 (0.31)
Daclizumab	–	1 (0.16)
Mycophenolate mofetil	–	2 (0.31)
Ocrelizumab	–	1 (0.16)
Rituximab	–	14 (2.2)
Missing	–	112 (17.58)

Initial analyses compared the distribution of within-person and between-person dissimilarity in overall metabolomic profiles; profiles were largely consistent within individuals; the distribution of dissimilarity between subjects was significantly larger than within-subject dissimilarity, as expected (p < 1E−16; Figure S1). The distributions of within-person and between-person Mahalanobis dissimilarity were similar when stratified by disease status (e.g., MS versus HC).

Metabolomic differences between MS and HC

Metabolic dysfunction score in overall metabolomic profiles

Overall, we detected shifts in age- and gender-adjusted metabolomic profiles in PwMS relative to HC (Figure 1A; p = 4.89E−7) using our age- and gender-adjusted metabolomic dysfunction classifier score. We identified samples that were highly divergent from the reference population (defined as >90th percentile of the metabolomic dysfunction score in HC); 131 (20.56%) MS patients had metabolomic dysfunction scores greater than this threshold, providing additional evidence of divergence in overall metabolite profiles (p = 7.65E−5). Differences in overall metabolomic dysfunction were consistent when we only included patients who were not on a disease-modifying therapy (DMT) at the time of blood collection (Figure S3).

Figure 1

Differences in metabolomic profiles between MS and healthy control (HC)

(A) Differences in overall metabolomic profiles. The dotted line denotes the 90th percentile of our metabolomic dysfunction classifier in HC.

(B) Volcano plot of results of individual metabolites. Each point denotes 1 metabolite. The x axis denotes the mean difference in metabolite levels between MS patients and HC, while the y axis denotes the −log(p value) for a test of the difference between MS and HC. Red-colored and labeled metabolites denote those significantly different between MS and HC (FDR-adjusted p < 0.05).

(C) The distribution of metabolites that significantly differed between MS versus HC (FDR-adjusted p < 0.05); plotted distributions are scaled to the median within HC.

(D) Largest connected component of the agnostic metabolite network derived using weighted gene correlation network analysis (WGCNA). Each node represents a metabolite; the color of the node represents the module color, as labeled by WGCNA. The brown, green, yellow, pink, and red modules are significantly different between MS patients and HC.

(E) The blue network represents all Metacyc metabolic reactions; nodes represent a metabolite and edges represent a metabolic reaction (e.g., reactant-product connection). We mapped results of the individual metabolites that differed between MS versus HC with FDR < 0.25 and extracted subnetworks of enriched metabolite interactions (corresponding roughly to metabolic pathways) using Prize-collecting Steiner Forest (PCSF) graph optimization. The bottom networks represent the 2 largest subnetworks extracted after PCSF optimization. All of the plots in this figure were derived from n = 984 samples.

See also Figures S2–S4 and Data S1B–S1D, S1F, and S1G.

Differences in individual metabolites and metabolic pathways

In analyses using individual metabolites, we found strong reductions among MS patients in lactate-related metabolites in aromatic amino acid (AAA) pathways. For example, phenyllactate (PLA), 3-(4-hydroxyphenyl)-lactate, indolelactate (ILA), and imidazole lactate were significantly reduced in PwMS (all false discovery rate [FDR]-adjusted p < 1E−8; Figures 1B and 1C; Table 2). Results were consistent when comparing HC to PwMS in analyses where we (1) included only MS patients who were not on a DMT at the time of blood collection, (2) included only progressive MS (PMS) patients, (3) excluded patients for whom treatment status was missing, and (4) restricted it to MS patients on any DMT. For example, PLA, 3-(4-hydroxyphenyl)-lactate, ILA, and imidazole lactate were reduced in PwMS versus HC regardless of the subset of patients included (Data S1B). Findings were also consistent in sensitivity analyses excluding pediatric samples. Similarly, in analyses in which we grouped metabolites a priori based on related biological functions (e.g., glutathione metabolism, tryptophan metabolism), we identified highly statistically significant differences in metabolic pathways including other AAA metabolites (Figures 1D and 1E; Table 3; all FDR-adjusted p < 1E−4). These results were consistent across different types of pathway and network analyses (Table 3; Data S1C). We also found strong reductions in branched-chain amino acid (BCAA)-related metabolites in both individual and pathway-based analyses. We observed differences in metabolite pathways related to bile acid metabolism (FDR-adjusted p = 3.14E−3; Table 3), xanthine metabolism (FDR-adjusted p = 4.91E−5), and acetylated amino acid metabolism (FDR-adjusted p = 4.19E−4). In analyses using compounds and reactions from the MetaCyc reaction database (as a proxy for estimated enzyme activity), we also identified a network of amino acids (predominantly AAA) and xanthine (caffeine) metabolites enriched in PwMS. Likewise, we also noted alterations in ratios of tyrosine (reactant) to 4-hydroxyphenylpyruvate (product) and phenylalanine to phenylpyruvate in PwMS, suggesting a reduced activity of enzymes involved in these metabolic reactions (Data S1D). Intriguingly, in analyses comparing age-matched PMS patients to relapsing remitting patients, differences in 3-(4-hydroxyphenyl)-lactate, ILA, phenylacetylglutamine, and p-cresol glucuronide were generally more extreme in PMS patients (Figure 2B; Data S1E).

Table 2

Analyses of metabolomic differences between multiple sclerosis (MS) and healthy controls (HCs) for individual metabolites with FDR-adjusted p < 0.05

Metabolite	Metabolite pathway	HMDB IDa	Metabolite modulea	Mean difference between MS versus HC (95% CI)b	p	FDR-adjusted p
Phenyllactate (PLA)	phenylalanine metabolism	HMDB00779	pink	−6.54 (−7.93 to −5.15)	3.16E−20	8.49E−18
3-(4-Hydroxyphenyl)lactate	tyrosine metabolism	HMDB00755	pink	−6.19 (−7.62 to −4.76)	1.73E−17	4.63E−15
Indolelactate	tryptophan metabolism	HMDB00671	pink	−5.76 (−7.19 to −4.33)	2.61E−15	6.97E−13
Tryptophan	tryptophan metabolism	HMDB00929	brown	−4.51 (−5.87 to −3.15)	7.21E−11	1.92E−8
Imidazole lactate	histidine metabolism	HMDB02320	pink	−4.82 (−6.28 to −3.35)	1.12E−10	2.98E−8
4-Hydroxyphenylpyruvate	tyrosine metabolism	HMDB00707	brown	−4.47 (−5.89 to −3.06)	5.89E−10	1.55E−7
α-Hydroxyisocaproate	leucine, isoleucine, and valine metabolism	HMDB00746	pink	−4.29 (−5.74 to −2.85)	6.22E−9	1.64E−6
Asparagine	alanine and aspartate metabolism	HMDB00168	brown	−4.42 (−5.93 to −2.90)	1.03E−8	2.70E−6
Methionine	methionine, cysteine, S-adenosyl methionine (SAM), and taurine metabolism	HMDB00696	brown	−3.97 (−5.36 to −2.59)	1.86E−8	4.86E−6
N-Acetylcarnosine	dipeptide derivative	HMDB12881	pink	−3.83 (−5.18 to −2.49)	2.36E−8	6.13E−6
2-Hydroxy-3-methylvalerate	leucine, isoleucine, and valine metabolism	HMDB00317	pink	−3.98 (−5.51 to −2.46)	3.03E−7	7.85E−5
Tyrosine	tyrosine metabolism	HMDB00158	brown	−3.51 (−4.89 to −2.12)	7.25E−7	1.87E−4
Choline phosphate	phospholipid metabolism	HMDB01565	gray	−3.42 (−4.80 to −2.04)	1.20E−6	3.07E−4
N-Acetylleucine	leucine, isoleucine, and valine metabolism	HMDB11756	green-yellow	−3.35 (−4.75 to −1.95)	2.60E−6	6.67E−4
Catechol sulfate	benzoate metabolism	HMDB59724	turquoise	−3.36 (−4.77 to −1.95)	2.98E−6	7.59E−4
Histidine	histidine metabolism	HMDB00177	brown	−3.15 (−4.51 to −1.80)	5.07E−6	1.29E−3
Phenylpyruvate	phenylalanine metabolism	HMDB00205	pink	−3.24 (−4.64 to −1.84)	5.71E−6	1.45E−3
N-(2-Furoyl)glycine	food component/plant	HMDB00439	turquoise	−3.38 (−4.87 to −1.89)	8.84E−6	2.23E−3
Isovalerylcarnitine (C5)	leucine, isoleucine, and valine metabolism	HMDB00688	red	−3.13 (−4.55 to −1.70)	1.71E−5	4.29E−3
Glucuronate	aminosugar metabolism	HMDB00127	gray	2.82 (1.53–4.11)	1.89E−5	4.73E−3
Theobromine	xanthine metabolism	HMDB02825	gray	−2.80 (−4.09 to −1.51)	2.06E−5	5.12E−3
Hippurate	benzoate metabolism	HMDB00714	turquoise	−3.06 (−4.47 to −1.65)	2.12E−5	5.26E−3
O-Methylcatechol sulfate	benzoate metabolism	HMDB60013	turquoise	−2.99 (−4.40 to −1.57)	3.60E−5	8.90E−3
Tiglylcarnitine (C5:1-DC)	leucine, isoleucine, and valine metabolism	HMDB02366	red	−3.12 (−4.61 to −1.64)	3.66E−5	9.00E−3
3-Hydroxybutyrate (BHBA)	ketone bodies	HMDB00357	blue	2.82 (1.47–4.16)	3.92E−5	9.60E−3
2-Methylbutyrylcarnitine (C5)	leucine, isoleucine and valine metabolism	HMDB00378	red	−3.03 (−4.49 to −1.58)	4.52E−5	1.10E−2
N-Palmitoylglycine	fatty acid metabolism (acyl glycine)	HMDB13034	blue	2.92 (1.48–4.37)	7.22E−5	1.75E−2
Betaine	glycine, serine, and threonine metabolism	HMDB00043	gray	−2.87 (−4.33 to −1.40)	1.22E−4	2.95E−2
p-Cresol-glucuronide∗	tyrosine metabolism	HMDB11686	purple	2.80 (1.36–4.23)	1.32E−4	3.18E−2
3-Methyl catechol sulfate (1)	benzoate metabolism	HMDB0240662	turquoise	−2.77 (−4.21 to −1.32)	1.72E−4	4.13E−2
N-Acetyltyrosine	tyrosine metabolism	HMDB00866	green-yellow	−2.73 (−4.15 to −1.30)	1.81E−4	4.33E−2
Xanthurenate	tryptophan metabolism	HMDB00881	pink	−2.82 (−4.30 to −1.34)	1.89E−4	4.51E−2

HMDB, Human Metabolite Database identifier. Metabolite module denotes the WGCNA module in which the metabolite falls.

Mean differences are adjusted for age, gender, and race; values displayed denote mean differences per SD in metabolite levels.

Table 3

Analyses of metabolomic differences between multiple sclerosis (MS) and HCs using WGCNA

Modulea		Module sizeb	Including all samples (n = 954)			Restricting to MS patients on no therapy and HCs (n = 434)			Restricting to MS patients on any therapy and HCs (n = 710)
	Overview of module composition		Mean difference between MS versus HC (95% CI)c	p	FDR-adjusted p	Mean difference between MS and HC (95% CI)c	p	FDR-adjusted p	Mean difference between MS and HC (95% CI)c	p	FDR-adjusted p
Pink	aromatic amino acids, branched chain amino acids	12	−2.12 (−2.58 to −1.65)	2.13E−19	2.77E−18	−1.87 (−2.51 to −1.22)	1.55E−8	2.02E−7	−1.95 (−2.48 to −1.43)	3.39E−13	4.41E−12
Brown	other amino acids	22	−1.16 (−1.64 to −0.69)	1.29E−6	1.55E−5	−1.43 (−2.09 to −0.78)	1.65E−5	1.99E−4	−0.79 (−1.32 to −0.26)	3.60E−3	3.24E−2
Turquoise	xanthine metabolites	26	−0.97 (−1.43 to −0.50)	4.17E−5	4.59E−4	−0.96 (−1.67 to −0.24)	8.75E−3	8.75E−2	−1.33 (−1.88 to −0.79)	1.86E−6	2.24E−5
Green-yellow	acetylated amino acid derivatives	7	−0.96 (−1.41 to −0.50)	4.24E−5	4.59E−4	−0.81 (−1.53 to −0.09)	2.65E−2	2.39E−1	−1.10 (−1.62 to −0.59)	2.94E−5	3.23E−4
Red	branched chain amino acids; carnitine metabolites	14	−0.93 (−1.39 to −0.47)	7.94E−5	7.15E−4	−0.96 (−1.63 to −0.28)	5.55E−3	6.10E−2	−0.91 (−1.43 to −0.39)	5.60E−4	5.60E−3
Magenta	bile acid metabolites	9	−0.78 (−1.25 to −0.32)	9.36E−4	7.48E−3	−0.75 (−1.45 to −0.05)	3.52E−2	2.81E−1	−0.66 (−1.19 to −0.13)	1.40E−2	9.79E−2
Blue	fatty acid metabolites; acyl carnitine metabolites	25	0.52 (0.05–0.99)	2.97E−2	2.08E−1	0.68 (−0.03 to 1.39)	5.87E−2	4.11E−1	0.27 (−0.25 to 0.80)	3.07E−1	4.49E−1
Tan	vitamin metabolites	5	−0.31 (−0.75 to 0.14)	1.81E−1	9.03E−1	−0.00 (−0.63 to 0.62)	9.94E−1	1.00E+0	−0.54 (−1.07 to 0.00)	5.08E−2	2.54E−1
Yellow	steroid metabolites	21	0.12 −0.34 to 0.58)	6.11E−1	1.00E+0	0.43 (−0.23 to 1.09)	2.01E−1	1.00E+0	0.61 (0.12–1.11)	1.51E−2	9.79E−2
Purple	aromatic amino acids, benzoate metabolism	7	−0.07 (−0.53 to 0.38)	7.48E−1	1.00E+0	0.07 (−0.62 to 0.75)	8.46E−1	1.00E+0	−0.43 (−0.96 to 0.10)	1.12E−1	4.49E−1
Black	gamma glutamyl metabolites	12	0.06 (−0.38 to 0.51)	7.80E−1	1.00E+0	−0.28 (−1.00 to 0.44)	4.48E−1	1.00E+0	0.39 (−0.11 to 0.89)	1.28E−1	4.49E−1
Green	methionine, glycine metabolism-related metabolites	16	−0.05 (−0.51 to 0.41)	8.30E−1	1.00E+0	−0.18 (−0.82 to 0.45)	5.76E−1	1.00E+0	0.38 (−0.12 to 0.88)	1.36E−1	4.49E−1

Module name denotes arbitrary name of set of metabolites identified by WGCNA.

Module size denotes the number of metabolites in a given module.

Mean differences are adjusted for age, gender, and race; values displayed denote mean differences in eigenmetabolite levels scaled by a factor of 102.

Figure 2

Association between metabolomic profiles and disability in people with MS

(A) Differences in overall metabolomic profiles by disability status (HC versus no cane versus use of cane). The dotted line denotes the 90th percentile of metabolomic dysfunction in HC.

(B) Volcano plot of results of comparison of individual metabolite concentrations between relapsing-remitting MS (RRMS) and age-matched progressive MS (PMS) participants; each dot is a single metabolite, and the x axis denotes the mean difference between RRMS and age-matched PMS participants, while the y axis denotes the −log(p value) for a test of difference between the 2 groups. Red-colored and labeled metabolites denote those significantly different between MS and HC (FDR-adjusted p < 0.05).

(C) Volcano plot of results of individual metabolites; the x axis denotes the association between metabolite level and EDSS scores, while the y axis denotes the −log(p value) for a test of association between metabolite level and EDSS score. The dotted horizontal line denotes nominal significance (e.g., −log[p = 0.05]).

(D) The distribution of metabolites that were potentially associated with EDSS scores by disability status (FDR-adjusted p < 0.15); plotted distributions are scaled to the median within HC.

(E) Oxidative and reductive pathway metabolism of aromatic amino acids (adapted from Dodd et al.). Metabolites colored in blue denote reductive pathways metabolites (and are the denominators for ratios plotted included in Figure 3F). Metabolites colored in red denote oxidative pathway metabolites (and are the numerators for ratios plotted in Figure 3F).

(F) The association between oxidative versus reductive pathway metabolism for selected aromatic amino acids (AAAs) metabolites. For each test, we modeled the ratio of oxidative terminal metabolites to reductive terminal metabolites, as identified from previous studies and noted in Figure 3D. The top plot for each ratio denotes differences between MS patients and HCs, while the bottom panel denotes the difference in the ratio of metabolites between MS patients using a cane (e.g., having a high level of disability) relative to those not using a cane (e.g., lower levels of disability). p values are derived from generalized estimation equation (GEE) models adjusting for age, gender, and race. For models of disability status, we evaluated the association between metabolite ratios and continuous EDSS using a similarly adjusted GEE model.

All of the plots in this figure were derived from n = 417 samples, except for (B), which was derived from n = 260 samples. See also Figure S4 and Data S1D, S1E, and S1H–S1L.

Figure 3

Differences in pathway activity scores between MS and HC by cell type

(A) Heatmap of amino acid pathway activity differences between MS and HC across cell types. The color of the cell denotes the degree of difference in pathway activity scores; darker blue cells denote pathway activity scores that are higher in HC relative to MS, whereas darker orange cells denote pathway scores that are higher in MS relative to HC. ∗p < 0.05, ∗∗p < 0.01.

(B) Differences in pathway activity scores for monocyte-like cell clusters for AAA pathways between MS and HC samples from CSF (left; n = 11 samples) and blood (right; n = 10 samples).

(C) Differences in network activity scores for monocyte cell clusters for AhR and HCA3 signaling pathways between MS and HC samples (n = 21 samples). ∗p < 0.05, ∗∗p < 0.01.

See also Figure S5 and Data S1M.

Sensitivity analyses assessing whether observed differences in metabolites and metabolic pathways in primary analyses were associated with different MS DMTs

For metabolites that were significantly different between MS patients and HC (e.g., those metabolites with FDR-adjusted p < 0.05 between MS patients and HC), we also assessed whether observed metabolomic differences were related to MS DMTs by comparing DMT classes (based on efficacy) versus no therapy and individual DMTs (among DMTs with ≥30 users) versus no therapy. For selected metabolites with significant differences between MS and HC, individual DMTs or DMT class did not appear to be strongly associated with differences in levels for these metabolites (Data S1F).

Other sensitivity analyses for primary metabolomics analyses between MS versus HC

The adjusted mean difference between MS and HC in a given metabolite level was not associated with how strongly a metabolite correlated with age in HC (correlation [βMS versus HC, βAge in HC] = −0.09; p = 0.16; Figure S4A). We also observed consistent results when stratified by median age (<40 years, ≥40 years; Data S1G). Results were relatively consistent when we adjusted for BMI when available. Results were also consistent when repeating analyses using a leave-one-out procedure excluding individual batches and repeating all of the analyses in the remaining set. Results were also similar in stratified analyses by serum versus plasma.

Association of metabolomics profiles with MS characteristics and disease severity measures

Having established that significant alterations exist in the circulating metabolome in MS patients compared to HCs, we then addressed the question of whether alterations in the metabolome are linked to MS disease severity (clinical or based on imaging [optical coherence tomography, OCT] measures). Higher levels of metabolic dysfunction were associated with increased disability status; MS patients using a cane had the highest dysfunction scores (Figure 2A). In analyses assessing the association between individual metabolites and EDSS scores, we observed that reductions in several AAA metabolites were associated with higher EDSS scores independent of age (Figures 2C and 2D; Data S1H). Results were consistent when restricting to subsets of MS patients in which we (1) included only patients who were not on a DMT at the time of blood collection, (2) included only PMS patients, (3) excluded patients for which treatment status was missing, and (4) restricted to MS patients on any DMT (Data S1H). Specifically, reduced levels of 3-(4-hydroxyphenyl)lactate (tyrosine metabolism) were significantly associated with higher EDSS scores (FDR-adjusted p < 0.05). Reductions in other AAA metabolites (ILA, imidazole lactate, PLA, kynurenine, kynurenate, tryptophan, and phenylalanine) were also nominally associated with higher EDSS scores (unadjusted p < 0.05; Data S1H). Several of these metabolites are largely produced by the reductive pathway of gut microbial AAA metabolism. In network or pathway-based analyses, we also observed consistent significant associations between AAA pathway metabolites and EDSS scores and in analyses stratified by age (Data S1I–S1K). In addition, increases in multiple gut microbiota-derived AAA metabolites (which are also uremic metabotoxins)—p-cresol glucuronide, p-cresol sulfate, and phenylacetylglutamine—were associated with greater disability (Figures 2C and 2D). These metabolites are derived from the oxidative pathway of gut microbial AAA metabolism. To assess whether the balance between the reductive and oxidative metabolism of AAA in the gut was associated with MS and disease severity, we used ratios of metabolites derived from these pathways for individual amino acids (Figures 2E and 2F). We noted highly significant associations between the ratio of oxidative to reductive metabolites of phenylalanine (phenylacetylglutamine to PLA ratio), tyrosine (p-cresol glucuronide or sulfate to 3-(4-hydroxyphenyl)lactate), and tryptophan (indole acetate to ILA ratio) and EDSS scores, suggesting that a shift in reductive versus oxidative metabolism of AAA may be related to disability (Figure 2F). In addition, results were consistent when we used the Age-Related Multiple Sclerosis Severity (ARMSS) (rather than the EDSS) as a measure of disease severity. Sensitivity analyses also confirmed that the association between a given metabolite and EDSS was not also associated with how strongly a metabolite correlated with age in HC (Figure S4B).

A subset of participants had OCT-derived measures of disease severity available. In these participants, AAA metabolites (xanthurenate, p-cresol-glucuronide∗, 4-hydroxyphenylpyruvate, phenylacetylglutamine, and PLA) and metabolite ratios (phenylacteylglutamine to PLA and p-cresol glucuronide to 3-(4-hydroxyphenyl)lactate) were also associated with differences in ganglion cell + inner plexiform layers (GCIPL) thickness (Data S1L). Thus, results were consistent using either clinical or imaging measures of disease severity.

Interestingly, we also noted a significant positive association between the magnitude of the difference between MS and HC for a given metabolite and the magnitude of the association between individual metabolite level and EDSS scores (correlation [βMS versus HC, βEDSS] = 0.41; p = 1.35E−12). For example, the metabolite p-cresol-glucuronide is elevated in PwMS relative to HC; higher levels are also associated with higher levels of disability (Figures 1B and 2C).

Single-cell metabolic gene expression in MS versus HC in blood and CSF

We used publicly available single-cell RNA sequencing (scRNA-seq) data from an existing study comparing cellular-level differences in metabolic gene expression from an integrative analysis of blood and CSF from MS patients and HC. After downloading raw cell counts, we implemented a standard quality control and cell clustering procedure to identify clusters of cells using marker gene expression (Figure S5) for use in subsequent analyses. Because of the strong differences in amino acid metabolism (specifically in AAA metabolism) identified using circulating metabolomics, we concentrated our scRNA-seq analyses of pathway differences in cell-type-specific gene expression between MS and HC in blood and CSF on these metabolic pathways in particular (Figure 3A). We identified significant differences in the pathway activity of AAA metabolic pathways between MS and HC in monocyte cell clusters (Figure 3B), including both clusters of cells enriched in blood and CSF; PwMS tended to have significantly lower levels of pathway activity in theses cell types. Relatedly, as we hypothesized that the observed abnormalities in AAA metabolism may also imply an altered balance of immunomodulatory metabolites (e.g., AhR or HCA3 agonists), we also tested whether a network of genes interacting with AhR or HCA3 was similarly altered; networks were defined using protein-protein interaction databases. We detected lower levels of AhR network genes in monocytes in PwMS relative to HC (p < 0.001; Figure 3C; Data S1M). We also found lower levels of expression for HCA3-network genes in PwMS; however, this difference was not statistically significant (p = 0.11; Figure 3C; Data S1M).

Effects of AAA-derived metabotoxins and lactate metabolites on human monocytes

Besides serving as biomarkers, circulating metabolites can have direct effects on cells through a variety of receptors. We first tested the effects of AAA-derived metabotoxins—indole acetate and phenylacetyl glutamine—that were identified as being related to higher EDSS severity and lower GCIPL thickness, on human peripheral blood mononuclear cells (PBMCs) from healthy controls (Figure 4A). Since our scRNA-seq analyses implicated potential differences in AAA metabolism in monocytes, we concentrated our analyses on this cell type. We noted an increase in tumor necrosis factor-α (TNF-α) production from CD14high monocytes with indole acetate (IAA) treatment with a clear dose-response relationship (Figures 4B and 4C). In a subset of participants, we also noted an increase in interleukin-6 (IL-6) production from CD14high monocytes with either phenylacetylglutamine (PAG) or IAA treatment (Figure S6B).

Figure 4

Effects of AAA-derived metabolites on human monocytes

(A)We treated freshly isolated peripheral blood mononuclear cells from HCs with either vehicle or escalating doses of AAA-derived metabotoxins (in duplicate), identified from our initial metabolomics analyses, for 24 h.

(B). We then used multiparametric flow cytometry to evaluate the change in pro-inflammatory cytokine production from monocytes (CD14high) and noted an increase in TNF-α production with increasing doses of IAA; a representative plot is shown here.

(C) Quantification of the increase in the proportion of TNF-α-producing monocytes with IAA treatment (n = 9). For all boxplots, the center line indicates the median, the box indicates the 25th and 75th percentiles, the whiskers indicate 1.5× interquartile range (IQR), and the dots indicate outliers.

(D) We isolated CD14+ monocytes from HCs and treated them with either vehicle or escalating doses of AAA-derived lactate metabolites (in duplicate), either overnight and then evaluated endocytosis of flurophore-labeled dextran, or for 36 h and then assessed cytokine production.

(E) We used flow cytometry to evaluate endocytosis of AF647-dextran by monocytes and noted an increase with indolelactate (ILA) treatment; a representative plot is shown here.

(F) Quantification of proportion of monocytes endocytosing AF647-dextran with ILA treatment (n = 5). The data for the various ILA doses is normalized to vehicle for each individual. Bars represent means and error bars represent standard errors of the mean.

(G) Boxplot of proportion of monocytes producing IL-6 with vehicle or various doses of ILA.

(H) Boxplot of proportion of monocytes producing IL-1β with vehicle or various doses of ILA. ∗p < 0.05, ∗∗p < 0.01, and ∗∗∗∗p < 0.001.

Also see Figure S6.

Since we noted that alterations in the balance of oxidative and reductive-pathway AAA metabolites are related to disease severity, we next assessed the effects of AAA-derived lactate metabolites ILA and PLA on human monocyte (CD14+) function (Figure 4D). We noted that treatment with ILA 50 μM led to an increased proportion of cells endocytosing dextran-AF647 (Figures 4E and 4F) compared to vehicle. We did not note an effect of PLA on monocyte endocytosis. We also noted that treatment with ILA 50 μM reduced IL-6 and IL-1β production from monocytes compared to vehicle (Figures 4G and 4H).

Discussion

This large-scale study identified differences in metabolomic profiles in PwMS relative to HC using a set of individual and pathway-based analyses. Several AAA were altered in MS and lower levels of AAA metabolites derived from the reductive metabolic pathway, and increased levels of AAA metabolites derived from the oxidative pathway were associated with higher disability scores. In addition, analysis of scRNA-seq data from blood- and CSF-derived immune cells demonstrated altered AAA metabolism in MS compared to controls. Finally, the identified AAA metabolites linked to disability status had functional effects on human monocytes modulating pro-inflammatory cytokine production and endocytosis in this population of immune cells.

We identified a marked disruption of multiple amino acid metabolic pathways, especially in the downstream metabolism of several AAAs that are primarily derived from gut microbial reactions (Figures 2E and 2F). Specifically, we noted a broad shift in AAA toward oxidative pathway metabolites relative to reductive pathway metabolites (e.g., Figure 2F). For example, reductive pathway metabolites (e.g., 3-[4-hydroxy]phenylactate [tyrosine], indole lactate [tryptophan], PLA [phenylalanine]) were largely reduced in PwMS. These specific AAA metabolites are found in large quantities in fermented foods and were recently identified as agonists for the HCA3 receptor, which is highly expressed on innate immune cells in humans and is thought to be involved in anti-inflammatory responses. Furthermore, these metabolites (e.g., ILA) can serve as agonists for the AhR and mediate immunosuppressive actions in immune and glial cells (e.g., limiting pathogenic activities of astrocytes); lower levels could result in increased inflammatory disease activity, both in the peripheral immune system and within the CNS. Consistent with observations that an imbalance between these two pathways has a role in MS, we noted strong associations between oxidative AAA pathway metabolites and increased MS risk and disease severity. Intriguingly, these results were consistent with analyses of metabolic pathway-level changes in gene expression in several cell types, most notably in a cluster of CSF-enriched monocytes whose marker gene signature resembled homeostatic microglia.,

Phenylacetylglutamine, derived from phenylacetate, is a product of the gut microbial oxidative metabolism of phenylalanine., A recent study demonstrated that phenylacetylglutamine (phenylalanine oxidative product) is associated with the risk of cardiovascular disease due to a direct effect on platelet reactivity through adrenergic receptor signaling. Similarly, the products of oxidative metabolism of tyrosine (p-cresol sulfate and p-cresol glucuronide) and tryptophan (indole acetate) are also metabotoxins associated with inflammation promotion. They are altered in other neurologic diseases and can increase the risk for cardiovascular disease and cardiac dysfunction. Several epidemiologic studies in MS have demonstrated a negative impact of vascular and related comorbidities on MS disease severity, and these metabolic alterations may serve as a link between the two conditions., Intriguingly, specific bacterial genes are required for reductive versus oxidative (porA for the oxidative pathway and fldH for the reductive pathway); future studies will explicitly link shifts in AAA with changes in these candidate metagenomic genes. Functional analyses also demonstrated that treatment with indole acetate led to increased production of the pro-inflammatory cytokine TNF-α from human monocytes. In additional functional analyses, we noted that ILA led to reduced production of the pro-inflammatory cytokines IL-6 and IL-1β and also promoted endocytosis (a function generally associated with an anti-inflammatory phenotype). Overall, this suggests that the imbalance in circulating levels of oxidative (higher) and reductive (lower) pathway AAA metabolites in MS would lead to more inflammatory myeloid cells producing larger amounts of key pro-inflammatory phenotypes. TNF-α, IL-6, and IL-1β have important roles in MS disease pathogenesis,23, 24, 25 and the effect of metabolites on the production of these cytokines from myeloid cells provides a possible mechanism by which alterations in the metabolome could lead to worsened MS disease activity and severity.

Also consistent with prior studies, we noted changes in other tryptophan pathway metabolites, with lower circulating levels of both circulating tryptophan and its endogenous metabolism-derived (e.g., kynurenine) metabolites in MS compared to HC., As described above, some tryptophan metabolites limit CNS inflammation via AhR-mediated mechanisms in both microglia and astrocytes. In the experimental autoimmune encephalitis (EAE) model of MS, microglial AhR deletion worsened EAE, increasing demyelination and CNS monocyte recruitment. Levels of related amino acid metabolites were also reduced in PwMS; the source of these changes remains unclear: diet, altered gut microbiota metabolism of these amino acids, or increased consumption by activated immune cells all may contribute. For example, mechanistic studies have shown that in mice lacking the programmed cell death protein-1 (PD-1) inhibitory checkpoint receptor on T cells, systemic decreases in tryptophan and tyrosine were due to increased uptake by activated T cells. Downstream reductions in circulating levels of these specific amino acids led to a deficiency in the neurotransmitters serotonin and dopamine in the brain and resulted in increased anxiety-like behaviors. Mood disorders such as depression and anxiety are very common in PwMS; it is possible that reduced levels of these metabolites could be a contributing cause of these common comorbidities in MS.

Several other metabolic pathways, including BCAA, bile acid metabolism, and xanthine metabolism, were altered in MS compared to HCs and have been previously reported in smaller studies of PwMS or in animal models. Circulating levels of BCAAs were reduced in MS. BCAAs are critical for the maintenance of regulatory T cells and their suppressive function in vivo. Hence, lower BCAA levels could impair Treg number and function and predispose to increased inflammatory T cell activity. We recently identified reduced circulating bile acid levels in PwMS and found that bile acid supplementation prevented the polarization of astrocytes and microglia to neurotoxic phenotypes both in vitro and in EAE. Here, we also noted changes in xenobiotic metabolism, including xanthine/caffeine metabolism in MS. Caffeine metabolism is altered in other neurodegenerative disorders, and intake of caffeine has been linked to the risk of developing MS, as well as disease severity; our results may suggest that altered caffeine metabolism may underlie the observed associations. In complementary analyses, mapping identified metabolites to a network created from known metabolic reactions, we identified two networks of metabolites that were enriched in MS compared to controls – these again included a network of multiple AAA and BCAAs and a network of caffeine metabolites. Thus, both agnostic and pathway-based analyses yielded complementary results.

Prior studies in PwMS note divergent metabolomic profiles, which is consistent with our findings,,31, 32, 33, 34, 35, 36 and several studies note alterations in AAA-related metabolites., Intriguingly, one study also found reduced levels of ILA, which was linked with lower levels of ILA-producing bacteria in the gut microbiome in PwMS. Distinct from our findings, other metabolomics studies have noted differences in various lipid subspecies (phosphatidylcholine, short-chain fatty acids)., Differences between results from these studies and the results presented here may be related to differences in the analytical platform used to assess metabolomic profiles or differences in the underlying study populations. We also only included metabolites that were consistently measured across populations in our study. It is worth noting that the vast majority of previous studies were of relatively limited sample sizes, and a key strength of our study is its large size; that is, we are able to detect relatively precise estimates of metabolomic differences between PwMS and HC.

There are several noteworthy strengths of this study, which include its large sample size, multiple sites, detailed analyses, and comprehensive assessment of circulating metabolites that were performed using the same metabolomics platform. We also implemented a stringent QC protocol evaluating metabolite stability within- and between-persons and over time. We complemented our primary analyses with extensive sensitivity analyses; we repeated analyses using a leave-one-out procedure to ensure that results were not driven by a single site or batch. We were also able to integrate data from other omics sources in novel multi-omics analyses to help guide the interpretation of our findings.

In summary, we have demonstrated strong differences in metabolomic profiles between PwMS and HC individuals. Namely, we note consistent alterations across each AAA toward increased oxidative relative to reductive pathway activity (implicating a specific subset of gut microbial alterations) as being associated with both MS risk and disability. These changes suggest a shift toward the increased production of metabotoxins, with a potential reduction in the production of immunomodulatory metabolites. We also identify altered AAA metabolic gene expression in MS monocytic cell populations in the blood and CSF and demonstrate the direct effects of AAA-derived metabotoxins and lactate metabolites on human monocytes. Ultimately, this novel metabolomics study has implications both for advancing cross-omics methods to understand the disease and for providing critical new insights into candidate pathologic mechanisms contributing to MS.

Limitations of the study

Despite the strengths of the study noted above, there are some limitations. The cross-sectional design limits certain conclusions; longitudinal studies are needed to evaluate how metabolomic changes can influence MS risk or disability changes over time. We included individuals who are largely prevalent DMT users. Therefore, we could not optimally assess the effects of initiation of different DMTs on the metabolome; several smaller studies have identified changes in circulating metabolites associated with different DMTs. However, for metabolites that strongly differed between PwMS and HC, results were consistent when we (1) included only MS patients who were not on a DMT at the time of blood collection, (2) included only PMS patients, (3) excluded patients for whom treatment status was missing, and (4) restricted to MS patients on any DMT. It will be important for future longitudinal studies to systemically characterize the association between DMT initiation and resultant change to the metabolome. We also lacked comprehensive body mass index (BMI) information on all of the participants. Still, results were consistent after adjusting for BMI in the subset of individuals where this information was available. We also did not have information on the time of last meal for all of the participants. A detailed assessment of diet or other potentially relevant comorbidities affecting metabolomic profiles was also not available. Many key metabolites differing between MS patients and HCs are derived from gut microbial reactions, and we lacked this information. Future studies will link metabolomic profiles with both composition and metagenomic features of the gut microbiota. Also, several studies of PwMS note differences in short-chain fatty acids. Because of the specific methods required to detect and quantify these biomarkers accurately, we were unable to assess their association with MS in our cohort. While intriguing, analyses that identified shifts in metabolic gene expression in AAA pathways within specific cell types are cross-sectional. Thus, they cannot identify whether such changes are a contributing cause or a result of the observed differences in the circulating metabolome. We also note that metabolic changes in several of the cell types have not been extensively characterized in PwMS or HC.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Pavan Bhargava, MD (pbharga2@jhmi.edu).

Materials availability

This study did not generate new unique reagents.

Experimental model and subject details

Metabolomics study population

Study participants were pooled from three sites: Johns Hopkins MS Center (JHU), the University of California at San Francisco (UCSF), and the Henry Ford Hospital MS Center or Accelerated Cure Project (ACP; samples were acquired by Henry Ford).

JHU study population

JHU study participants (n = 640) were pooled from several ongoing or recently completed clinical research studies in which blood samples were acquired for metabolomics analyses. We included baseline metabolomics samples (e.g., before the initiation of the intervention) from two recently completed studies of dietary interventions in people with relapsing-remitting MS: one study included relapsing-remitting MS patients aged 18-75, with BMI ≥ 25 kg/m2 and receiving natalizumab, the other included relapsing-remitting MS patients, aged 18-50, untreated or treated with injectable MS therapies (interferon beta or glatiramer acetate) and a BMI ≥ 22.5 kg/m,, We also included samples from an observational study of relapsing-remitting MS patients who were treated with dimethyl fumarate and age- and race-matched HC. Finally, we included baseline samples from a study of Caucasian age-matched relapsing-remitting MS patients and HC who had serum 25-hydroxyvitamin D levels < 20ng/mL. For remaining studies, participants with MS (including both relapsing-remitting and progressive subtypes) and HC were recruited by convenience sampling from the JHU MS Center and provided blood samples. Participants for this study (both MS and HC) were members of a long-standing observational cohort and could contribute multiple samples; the MS samples were collected at different time points of treatment and were generally not pre- and post- treatment samples. Blood was processed within 3 h of collection using a standardized protocol that was common to all the studies conducted at JHU, and serum or plasma was aliquoted and stored at −80 C° until metabolomics analyses.

In a subset of study participants (n = 305), information on measures of disease severity, including measures of disability status (via the Expanded Disability Status Scale [EDSS]), were available. A subset of participants (n = 192) also had optical coherence tomography (OCT) scans, as an imaging biomarker of retinal neurodegeneration that were obtained within 6 months of metabolomics assessment. Segmentation of retinal layer thicknesses was performed using a validated automated segmentation algorithm, as previously described.57, 58, 59 We assessed the association between metabolite levels and thickness of the macular ganglion cell + inner plexiform layers (GCIPL); we selected GCIPL thickness as the primary OCT phenotype as this composite measure is more strongly associated with brain atrophy and is less vulnerable to swelling that may occur in the context of inflammation of the optic nerve., Eligible participants were those who did not have a history of diabetes mellitus, uncontrolled hypertension, glaucoma, prior ocular surgery or trauma, refractive errors exceeding ± 6 diopters, or other significant neurological or ophthalmological conditions.

Henry Ford/ACP study population

Individuals with MS (including both relapsing-remitting and progressive subtypes) as well as HC were recruited by convenience sampling from the Henry Ford Hospital MS Center. A subset of MS serum samples was acquired from serum repository at the Accelerated Cure Project (ACP). For both sets of MS samples, participants (n = 184 total) had to have a confirmed diagnosis of MS and provide blood samples, which were used for metabolomics assessment. Metabolomics assessment was conducted in a single batch, regardless of the original parent study (e.g., Henry Ford or ACP). Blood samples from HFH MS Center and ACP were obtained through an IRB approved study of the metabolomics signature in MS.

UCSF pediatric study population

MS patients and controls were selected randomly from two cohort studies of pediatric MS led by investigators at UCSF (n = 136). These cohorts included: 1) pediatric MS or clinically isolated (CIS) patients enrolled in a prospective cohort study from the UCSF Pediatric MS clinic, 2) pediatric MS or CIS patients who participated in a multi-center case-control study evaluating risk factors for pediatric MS and who also provided demographic and clinical information.,, For both cohorts, eligible MS participants were those with relapsing-remitting MS or CIS with high risk of conversion to MS (e.g., ≥ 2 T2 hyperintense foci on MRI), with disease onset < 18 years and who had seen a neurologist within 4 years of symptom onset. Diagnoses were confirmed by a pediatric MS specialist. HC were recruited from general or specialty pediatric clinics and were aged < 22 years, did not have a history of autoimmune disorders (except asthma or eczema) or previous history of severe health condition or treatment with immunosuppressive medications or have a biological parent or sibling with MS. All participants provided a blood sample at the time of recruitment into the study. Disability status via the EDSS score was also available on people with MS. Blood was processed and frozen at −80 C° within 3 h of collection.

Peripheral blood mononuclear cells from HC

We performed phlebotomy and collected blood from healthy controls (n = 9) and peripheral blood mononuclear cells (PBMCs) were isolated using SepMate tubes (STEM cell). These cells were then utilized for in vitro studies examining the effects of AAA-derived metabloites on human monocytes as detailed below.

Method details

Assessment of metabolomic profiles

All metabolomics analyses via mass spectrometry were conducted at Metabolon (Durham, NC), and methods have been described in detail elsewhere. JHU samples were pooled from metabolomics analyses conducted in 7 different batches, while UCSF and Henry Ford samples were each analyzed in separate single batches, totaling 9 batches overall. Briefly, samples were thawed and underwent additional preparation (derivatization), as previously described. The derivatized samples were subjected to either gas chromatography followed by mass spectrometry (GC/MS) or liquid chromatography followed by tandem mass spectrometry (LC/MS/MS). Mass spectra obtained from these techniques were then matched to a library of spectra derived from standards to identify specific metabolites, and the area under the curve for the mass spectra was used to calculate the relative abundance of each metabolite.

scRNA-seq of blood and CSF in MS versus HC

In follow-up analyses, we also assessed whether any metabolic differences we observed in the serum or plasma correlated with changes in gene expression at the cellular level. To do so, we used publicly available single cell RNA sequencing (scRNA-seq) data from an existing study comparing cellular-level differences in metabolic gene expression from an integrative analysis of blood and CSF from MS patients and HC. Eligible MS patients for the original study were treatment-naive patients with an initial episode suggestive of MS (CIS) or those with relapsing-remitting MS. Eligible control subjects were those providing blood and CSF samples following workup for idiopathic intracranial hypertension (IIH) matched to MS subjects by age, gender, and CSF features (protein, lactate, and glucose levels); full details are provided in Schafflick et al.

AAA-derived metabotoxins on human monocytes

To determine whether AAA-derived metabotoxins affect human immune cell function, we isolated peripheral blood mononuclear cells (PBMCs) from 9 healthy controls and cultured 5e5 cells in 96 well U-bottom plate for 24 h in AIM-V serum-free media (ThermoFisher) with either vehicle or varying doses of indole acetic acid (IAA) or phenylacetylglutamine (PAG) in duplicate. For the final 5 h of culture, monensin was added (BioLegend). At the end of this time period we washed cells then Fc blocked (TruStain, BioLegend) and stained with a fixable viability marker (Zombie NIR, BioLegend). We then stained with surface antibodies (CD14, CD16, CD83, CD86. BioLegend) followed by washing. Cells were then fixed and permeabilized overnight with IC Fix/Perm (ThermoFisher). Cells were washed in permeabilization buffer (ThermoFisher) and stained with cytokine markers (IL6, TNFα. BioLegend). Cells were then acquired on a flow cytometer (Cytek Aurora, Cytek Biosciences).

AAA-derived lactate metabolites on human monocytes

To determine whether AAA-derived lactate metabolites affect human monocytes, we isolated PBMCs from 5 healthy controls and isolated CD14+ monocytes by positive bead selection (CD14 UltraPure Microbeads, human, Miltenyi Biotec). CD14+ purity was typically ∼95% following isolation. Monocytes were then cultured in 200 μl serum-free defined media (AIM-V, ThermoFisher) at concentration of ∼2e5 cells/well. Cells were cultured either overnight (endocytosis assay) or for ∼36 h (cytokine assay) at 37°C, 5% CO2 in the presence of AAA-derived lactate metabolites (DL-Indole-3-lactic acid [ILA], Phenyllactic acid [PLA], Sigma-Aldrich) at specified concentrations or vehicle in duplicate. For endocytosis, after an overnight incubation 10kD dextran beads conjugated to Alexa Fluor 647 (ThermoFisher) were added at concentration of 50 μg/mL and incubated for 90 min at 4°C (negative control) or 37°C. Cells were then washed twice in FACS buffer (2% FBS, 2mM EDTA in PBS), stained with propidium iodide, and analyzed on a flow cytometer (Cytek Aurora, Cytek Biosciences). After gating on viable singlet cells, Alexa Fluor 647 fluorescence was evaluated. For cytokine assessment, monensin (Biolegend) was added for the last 5 h of the 36 h incubation. Cells were then Fc blocked and stained with a fixable viability marker (Zombie NIR, Biolegend) for 20 min at room temperature, then fixed with IC Fixation Buffer (ThermoFisher) for 30 min at room temperature and permeablized in Permeabilization Buffer (ThermoFisher). Intracellular staining was then performed in the permeabilization buffer for 1 h at room temperature. Cells were then washed twice with permeabilization buffer, once with FACS buffer, resuspended in FACS buffer and analyzed on a flow cytometer (Cytek Aurora, Cytek Biosciences).

Quantification and statistical analysis

Quality control

We initially included 329 metabolites that were measured in all 9 metabolomics runs in 960 total samples. We then implemented the following quality control (QC) procedure to identify potential outlying metabolites and/or samples. First, we removed metabolites with > 20% missing values across samples (n = 27) and imputed missing metabolite values using k-nearest neighbors (10 neighbors used for each imputation); consistent results were observed when using the minimum value of observed metabolites. We then log-transformed all metabolites. We adjusted for batch (and site/specimen type) using the ComBat algorithm, which is a harmonization technique that was designed to remove batch-related extraneous variation while conserving biological variation,; ComBat has been extended beyond gene expression (where it was originally developed) to adjust for scanner or site effects from segmented MRI volumetric data. Results before and after the application of ComBat by batch are displayed in Figure S2A and by sample matrix type - serum versus plasma in Figure S2B. Our metabolomics analyses also included both within-batch replicates (e.g., the same sample was included twice in the same metabolomics run; n = 70 samples) and between batch samples in which one serum and one plasma sample were included. We calculated the intraclass correlation coefficient (ICC) for both within- and between-batch replicate samples. ICCs for both within- and between-batch were generally high; the median ICC for within-batch replicates was 0.94 (IQR: 0.86, 0.98; Figure S2D) and the median for between-batch replicates was 0.79 (IQR: 0.63, 0.89; Figure S2C). We removed metabolites which had ICCs < 0.40 in 1) between batch replicates (n = 27 metabolites), 2) within-batch replicates (n = 4 metabolites) or 3) both (n = 2 metabolites). The final analyses included 269 metabolites. From metabolomic samples with within- or between-batch replicates, we randomly selected one from each set to be eligible for inclusion in the main analyses. We tested for potential sample outliers using principal components analyses (PCA); we excluded those who were > 3 standard deviations (SD) of the mean for PC1 and PC2 (n = 1 excluded). We also identified potential outliers using a Euclidean distance-based sample network (and a standardized network connectivity measure [Zk], as described in Langfelder and Horvath.); outlying samples were classified as those who had Zk values < −4 (n = 5 samples), as suggested by Langfelder and Horvath; leaving 954 samples (99.1%; n = 756 unique individuals) eligible for the analysis.

We also compared the distribution of within-person and between-person dissimilarity in overall metabolomic profiles, as several participants in our study contributed multiple samples. We calculated the within-person dissimilarity as the median Mahalanobis dissimilarity between metabolomic profiles taken from the same individual. To calculate between-person dissimilarity, we first averaged metabolomic profiles taken from the same individual to estimate a composite “average” metabolomic profile. We then calculated the median Mahalanobis dissimilarity between metabolomic profiles from different individuals.

Primary analyses comparing MS versus HC

For our primary analyses, we compared metabolomic profiles between people with MS and HC considering 1) global differences in the overall metabolome, 2) differences in individual metabolites, and 3) differences in composite metabolic pathway measures.

Metabolomic dysfunction classifier

Our initial goal was to test for overall differences in the circulating metabolome between people with MS and HC using a single summary measure which we denoted as a “metabolomic dysfunction” classifier to be used in statistical tests (as we hypothesized MS patients to have highly divergent metabolomic profiles and this measure allowed us to formally assess this hypothesis). To create this classifier, we set the ‘reference population’ as all age- and gender-adjusted metabolomic profiles from HC and derived a metabolomic dysfunction score as the median Mahalanobis dissimilarity among age- and gender-adjusted metabolomic profiles in MS patients to this reference set. To identify samples which were highly divergent from the reference population, we thresholded the metabolomic dysfunction score at the 90th percentile (e.g., samples which have metabolomic feature configurations that have < 10% probability of occurring in a person without MS. Sensitivity analyses also adjusted for race/ethnicity in metabolomic dysfunction scores.

Individual and pathway-based analyses

We compared individual log-transformed metabolites between people with MS and HC using generalized estimating equations (GEE; to account for multiple metabolomic profiles contributed by some participants) with a Gaussian link function and adjusted for age, gender, and race/ethnicity. Additional analyses adjusted for body mass index (BMI; kg/m2), which was only available for a subset of individuals (n = 568). For metabolic pathway-based analyses, we performed two sets of complementary analyses: an agnostic approach to discover sets of related metabolites and another incorporating known information on metabolic pathways and metabolite interactions. For the agnostic approach, we derived novel sets of related metabolites using a weighted gene-expression correlation network analysis (WGCNA); WGCNA is a systems biology-based technique that was originally developed to study correlation patterns in gene expression and has been extended to other settings (including our previous work in metabolomics), cancer, and neuroimaging.,, The goal of WGCNA is to find clusters (or metabolic modules in this case) of highly interconnected nodes (metabolites, here) using a correlation network. The identified metabolic modules can be summarized into single measures (i.e., eigen-metabolites) as the first principal component of the identified module and can be used in subsequent analyses. For a priori-based pathway analyses, we performed two sets of related analyses incorporating known biological information. For the first, we classified metabolites into groups (> 3 metabolites) based on related biologic function (e.g., glutathione metabolism, tryptophan metabolism, among others); pathway memberships for each metabolite are included as a part of Data S1B. We applied a resampling-based permutation-based algorithm to assign statistical significance while preserving metabolite-metabolite correlation. To do so, we fit individual models for each metabolite (via GEE, adjusting for age, gender and race/ethnicity, as above) and extracted and ranked the obtained p value; ranks for a given pathway set were then averaged. We then permuted phenotype labels 10,000 times and repeated the above procedure to calculate an average rank of the p values for a pathway set for each permutation. Finally, the p value for a given pathway set is the probability that the observed average rank is less than the expected average rank (as calculated from the 10,000 permutations). For the second, we used known metabolic reactions and interactions from available databases to create a metabolic network to use for subsequent analyses. To do so, we downloaded all MetaCyc compounds and corresponding reactions. Individual compounds served as the nodes, and we set an edge to be a known (or predicted) reactant-product metabolic interaction; in total, the network included 2248 nodes (metabolites) and 7449 edges (reactant-product metabolic interaction). We then mapped results of the individual metabolites that are potentially different between MS patients and HC (e.g., individual metabolites which differed between MS versus HC with FDR < 0.25) and extracted subnetworks of enriched metabolite interactions (corresponding roughly to metabolic pathways) using Prize-collecting Steiner Forest (PCSF) graph optimization. We also performed analyses assessing pairwise ratios of measured metabolites on opposite sides of MetaCyc-defined metabolic reactions, as a proxy for estimated enzyme activity of metabolic reactions. For this analysis, we included 139 metabolite-metabolite ratios participating in 207 metabolic reactions and tested for differences between MS and HC using multivariable-adjusted GEE, as above.

Sensitivity analyses

We performed an additional set of sensitivity analyses where repeated analyses using a leave-one-out procedure where we excluded individual batches and repeated all analyses in the remaining set to ensure that a single subgroup of samples was not driving the observed associations. Sensitivity analyses additionally stratified by serum and plasma origin of each individual sample and pooled results using meta-analysis methods.

Sensitivity analyses related to MS therapies

We also assessed whether metabolomic differences that we observed in primary analyses were associated with specific MS disease modifying therapies (DMTs) by comparing DMT efficacy classes (higher, moderate, lower) versus no therapy and individual DMTs (among DMTs with ≥ 30 users) versus no therapy. Higher efficacy medications included rituximab, natalizumab, daclizumab and mycophenolate mofetil. Moderate efficacy DMTs included dimethyl fumarate and teriflunomide, and lower efficacy DMTs included glatiramer acetate and interferon beta preparations.

Analyses of disease severity measures

We assessed the association between EDSS (417 EDSS scores from 312 individuals with MS) and metabolomic profiles (i.e., global metabolomic differences, individual metabolites, and pathway-based analyses) and metabolite ratios using multivariable-adjusted GEE, as above. In sensitivity analyses, we also assessed the age-related MS severity score (ARMSS; as a marker of MS severity) and metabolomic profiles (and did not adjust for age in these analyses). ARMSS is similar to the MS severity score, except it uses age instead of disease duration, as age is typically unbiased and more readily obtained than disease duration. In the subset of individuals with OCT imaging linked to a metabolomics sample (n = 192), we also assessed whether GCIPL thickness (as a quantitative imaging measure of retinal neurodegeneration) was associated with selected metabolites (or metabolite ratios) that demonstrated significant differences between MS and HC (metabolites with FDR-adjusted p < 0.05 among the 269 individual metabolite tests). Similar to above, we used GEE (to account for multiple eyes per person) and adjusted for a similar set of covariates. Similar sensitivity analyses additionally stratified by serum and plasma origin of each individual sample and pooled results using meta-analysis methods.

Pathway activity scores for scRNA-seq data

After downloading raw cell counts, we implemented a standard quality control procedure (e.g., excluding low quality samples with few detected cells and features, cells with fewer than 200 features [gene], cells with ≥ 2500 features, cells with < 5% of percentage of all features coming from mitochondrial gene sets,). Normalized samples were integrated using the Seurat pipeline v3 (e.g., selection of integration anchors using canonical correlation analysis and 30 PCs). Samples were then clustered using a shared nearest neighbor (SNN) modularity optimization-based clustering algorithm; clustering and sample integration were then inspected visually for quality and separation. Based on marker gene expression, we identified a several clusters of T cells: effector-memory (EM)-like CD4+ T cells (CD69), central-memory (CM)-like CD4+ T cells (CD27), naive CD4+ T cells (TRAC, CD4), activated CD8+ cells (CD8B and CCL5), non-activated CD8+ cells (CD8B and CCR7) and T-regulatory cells (FOXP3). We also identified clusters of NK cells (GNLY, NKG7), B cells (CD79A), and plasmablasts (IGHG). We also identified monocyte clusters and monocyte cell clusters that were mostly blood-derived (Mono-B; FCGR3A/CD16) or CSF-derived (Mono-CSF; CD14, C1QA, C1QB). Other myeloid lineage cells were clustered into mDC type 1 (mDC1; WDFY4, BATF3), mDC type 2 (mDC2; FCER1A, CD1C). Other clusters included plasmacytoid dendritic cells (pDC; TNFRSF21) and megakaryocytes (MegaK; GNG11). We detected one cluster of red blood cells (HBA1, HBA2, and HBB), which was removed from the analysis. tSNE plots following sample clustering and labeling are depicted for all cells in Figure S5A and by sample type (blood, CSF) in Figure S5B. Our primary goal for this analysis was to identify cell-type specific shifts in metabolic gene expression occurring between MS and HC at the pathway level. To do so, we calculated a cell-type specific pathway activity score using curated gene/pathways lists ways (including pathways with at least 10 genes) that are available in Xiao et al. and adapting a similar methodology described therein. Briefly, the pathway activity score is calculated through the following steps 1) calculate the mean expression level for gene for cell type for individuals as and the average gene expression across all cell types, 2) calculate average relative expression for a given cell type relative to the average expression across cell types as , 3) derive final pathway activity score for pathways as the weighted average of cell type-specific relative expression values, , for each participating gene as . Weights are derived as the number of pathways in which a gene participates. Pathway activity scores for each cell type were averaged for MS and HC and calculated overall and stratified by sample type (blood, CSF) separately. In analyses, we compared the ratio between MS and HC . Statistical significance was assessed by a permutation test of sample labels. As our primary results suggest that alterations in AAA metabolism may also imply an altered balance of immunomodulatory metabolites (e.g., AhR or HCA3 agonists that have been shown by previous studies to have relevance to MS), we also tested whether the network of genes interacting with AhR or HCA3 were similarly reduced in people with MS in specific cell types. To perform this test, we identified networks of genes interacting with AhR or HCA3 using the STRING protein-protein interaction (PPI) database, as it provides both physical and functional interactions (e.g., signaling cascades). We then tested for differences in network activity between people with MS and HC across cell types using a score derived similarly to the pathway activity scores. For example, the mean expression for gene in the network for each cell type is calculated. Then, the relative average expression of a gene in the network in a given cell type relative to its expression across all cell types is calculated. Finally, the network activity is determined by calculating a weighted average of cell type-specific relative expression values for each gene. The weights are determined by the number of edges in the overall STRING PPI. Similar to the pathway activity score, we calculated each network activity score (for AhR and HCA3 networks) in MS patients and in HC separately and compared the ratio between MS versus HC. Statistical significance was also determined via permutation.

We identified CD14high monocytes (gating strategy is shown in Figure S6A) – and determined the proportion producing pro-inflammatory cytokine production (IL6, TNF-a). We then analyzed the effects of these metabolites on human monocytes using a one-way ANOVA with Dunnett’s multiple comparisons test to compare proportion of cytokine producing cells between the different treatment groups.

For endocytosis assay we identified monocytes that were positive for Dextran-AF647 and measured the proportion of cells involved in endocytosis. We used a one-way ANOVA with Dunnett’s multiple comparisons test to compare proportion of monocytes endocytosing Dextran-AF647 between the different treatment groups. Similarly, we identified monocytes producing pro-inflammatory cytokines IL-6 and IL-1β and then utilized a one-way ANOVA with Dunnett’s multiple comparisons test to compare proportion of monocytes producing these cytokines between the different treatment groups.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Antibodies
CD14-FITC	Biolegend	Cat#367116; RRID: AB_2571929
CD16-PerCP	Biolegend	Cat#302030; RRID: AB_940380
CD83-PE	Biolegend	Cat#305308; RRID: AB_314516
CD86-BV711	Biolegend	Cat#305440; RRID: AB_2565835
TNFa-BV421	Biolegend	Cat#502932; RRID: AB_10960738
IL6-PECy7	Biolegend	Cat#501120; RRID: AB_2572042
IL1b-AlexaFluor 647	Biolegend	Cat#508208; RRID: AB_604135
Chemicals, peptides, and recombinant proteins
Indole-3-acetic-2,2-d2 acid	Sigma-Aldrich	Cat#492817-100MG; CAS#24420-86-8
Phenylacetyl L-Glutamine	Cayman Chemical Company	Cat#16724; CAS#28047-15-6
DL-3-Phenyllactic acid	Sigma-Aldrich	Cat#P7251-10G; CAS#828-01-3
DL-Indole-3-lactic acid	Sigma-Aldrich	Cat#I5508-1G-A; CAS#832-97-3
Deposited data
Metabolomics data deposited in Mendeley Data	Metabolon, Inc (Durham, NC)	Medeley Data: https://data.mendeley.com/datasets/zgtn6k2xsh/1
Single-cell RNA-seq data	GEO	GEO: GSE138266
Software and algorithms
Code to reproduce the main findings	Github	https://github.com/kfitzg13/metabolomics
R version 3.6.1	N/A	https://www.R-project.org/
sva version 3.34.0	N/A	https://bioconductor.org/packages/release/bioc/html/sva.html
ms.sev version 1.0.4	N/A	https://cran.r-project.org/web/packages/ms.sev/index.html
gdata version 2.18.0.	N/A	https://cran.r-project.org/web/packages/gdata/index.html
lubridate version 2.18.0	Grolemund and Wickham39	https://cran.r-project.org/web/packages/lubridate/index.html
foreign version 0.8-72	N/A	https://cran.r-project.org/web/packages/foreign/index.html
haven version 2.1.1	N/A	https://cran.r-project.org/web/packages/haven/index.html
impute version 1.60.0.	N/A	https://cran.r-project.org/web/packages/impute/index.html
ggplot2 version 3.3.3	Wickham40	https://cran.r-project.org/web/packages/ggplot2/index.html
WGCNA version 1.68	Langfelder and Horvath41,42	https://https://cran.r-project.org/web/packages/WGCNA/index.html
gee version 4.13-19	N/A	https://cran.r-project.org/web/packages/gee/index.html
geepack version 1.2-1	Højsgaard et al.,43 Yan and Fine,44 Yan45	https://cran.r-project.org/web/packages/geepack/index.html
scater version 1.14.6	McCarthy et al.46	https://cran.r-project.org/web/packages/scater/index.html
edgeR version 3.28.1	Robinson et al.,47 McCarthy et al.48	https://bioconductor.org/packages/release/bioc/html/edgeR.html
scran 1.14.6	Lun et al.49	https://bioconductor.org/packages/release/bioc/html/scran.html
uwot version 0.1.5	N/A	https://cran.r-project.org/web/packages/uwot/index.html
Rtsne version 0.15	N/A	https://github.com/jkrijthe/Rtsne
Seurat version 3.1.3	Hao et al.50	https://cran.r-project.org/web/packages/seurat/index.html
EnsDb.Hsapiens.v86 version 2.99.0	N/A	http://bioconductor.org/packages/release/data/annotation/html/EnsDb.Hsapiens.v86.html
BiocSingular version 1.2.2	N/A	https://github.com/LTLA/BiocSingular
BiocNeighbors version 1.4.2	N/A	https://bioconductor.org/packages/release/bioc/html/BiocNeighbors.html
SingleR version 1.0.5	Aran et al.51	https://bioconductor.org/packages/release/bioc/html/SingleR.html
pheatmap version 1.0.12	N/A	https://cran.r-project.org/web/packages/pheatmap/index.html
batchelor 1.2.4	Haghverdi et al.52	https://bioconductor.org/packages/release/bioc/html/batchelor.html
Matrix 1.2-17	N/A	https://cran.r-project.org/web/packages/Matrix/index.html
AUCell version 1.8.0	Aibar et al.53	https://bioconductor.org/packages/release/bioc/html/AUCell.html
	N/A
GSA version 1.03.1	N/A	https://cran.r-project.org/web/packages/GSA/index.html
ggrepel version 0.8.1	N/A	https://cran.r-project.org/web/packages/GSA/index.html
Flowjo	BD	10.7
Other
Dextran, Alexa Fluor 647, 10,000 MW, Anionic, Fixable	ThermoFisher	D22914
AIM V Medium, liquid (research grade)	ThermoFisher	12055091
CD14 MicroBeads UltraPure, human	Miltenyi Biotec	130-118-906