Metabolic alterations in multiple sclerosis and the impact of vitamin D supplementation

Multiple sclerosis (MS) is a chronic demyelinating disorder of the CNS that has both inflammatory and neurodegenerative components (1). Various “omics” technologies have improved our understanding of the pathogenesis of MS (2). Metabolomics refers to the measurement of multiple small molecules in biological media (e.g., plasma, serum, or cerebrospinal fluid) (3). Metabolomics provides information downstream of other techniques (e.g., genomics, transcriptomics) as well as insights into the gut microbiota and environmental exposures (4, 5). The few studies that have used metabolomics in MS either examined limited metabolites or were technically unable to definitively identify many metabolites (6–8).

Vitamin D insufficiency is a risk factor for developing MS and is linked to increased disease activity in those with established disease (9, 10). We recently demonstrated that MS patients have a lower rise in serum 25-hydroxyvitamin D [25(OH)D] levels compared with healthy controls (HCs), when given the same amount of oral cholecalciferol supplementation (11).

In other populations, vitamin D supplementation reduces markers of oxidative stress and also positively affects many other metabolic markers (12, 13). Whether MS patients display a similar change in oxidative stress after vitamin D supplementation is not known. Since oxidative stress is thought to contribute to neurodegeneration and, thus, long-term disability in MS, we sought to evaluate if global metabolomic differences, especially in metabolites related to oxidative stress, exist between MS patients and healthy individuals and whether changes in these markers could provide further biologic insight into the impact of 25(OH)D levels in people with MS.

Cross-sectional cohort

Metabolomic profiles differ between MS patients and HCs. The demographic and disease characteristics of the participants in the cross-sectional cohort (27 relapsing-remitting MS [RRMS], 27 HCs) are provided in Table 1. The overall design of the study is depicted in Figure 1. Metabolomic profiling detected 536 compounds, of which 492 (92%) passed quality control procedures and were eligible for analysis. After adjusting for age and sex, 43 metabolites differed between MS patients and HCs (P value cutoff < 0.05), with an average absolute standardized difference of 0.62 ± 0.20 in their relative abundance between groups. Using a partial least squares discriminant analysis (PLS-DA) model, the two groups were clearly separated as indicated by R²Y (a measure of goodness-of-fit of the model assessing the degree to which metabolic profiles differentiate MS versus HC) and Q²Y (a quality assessment measure estimating consistency between predicted versus original data estimated from cross-validation) statistics (Figure 2; R²Y = 0.95; [P = 0.047]; Q²Y=0.56 [P = 0.001]). Important metabolites driving separation between the groups were identified using high variable importance in projection (VIP) scores (that estimate the relative importance of a variable in the PLS model) and included a large set of metabolites involved in redox homeostasis (glutathione metabolism and γ-glutamyl amino acid metabolism) (Table 2). Other metabolites with high VIP scores were related to the urea cycle, histidine, and xenobiotic (benzoate, or caffeine) metabolism (Table 2).

Figure 1

Flow chart of study procedures. This figure depicts design of the study and steps beginning with sample collection for the two cohorts, followed by metabolomics analysis, preprocessing of metabolomics data, and then the various statistical analyses employed: partial least squares discriminant analysis, pathway analysis, and weighted correlation network analysis.

Figure 2

Cross-sectional cohort — partial least squares discriminant analysis of metabolomics data. (A) The plot depicts separation of the relapsing-remitting multiple sclerosis (RRMS) and healthy control (HC) groups (n = 27 in each group) utilizing the first two components of the partial least squares discriminant analysis (PLS-DA) model. (B) Overview of the PLS-DA model with four components; R²Y was 0.979 and Q²Y was 0.558. (C) Results of the 1,000-fold cross-validation of the PLS-DA model using permutation of the group labels to generate a distribution to derive P values. The P value for R²Y was 0.047 and for Q²Y was 0.001.

Table 1

Demographic characteristics of study populations

Table 2

Cross-sectional cohort VIP scores for metabolites driving separation between MS and HC groups in the PLS-DA model^A

Pathway analyses using WGCNA and a priori–defined pathway sets. Results of the weighted correlation network analysis (WGCNA) analysis classified metabolites into 15 modules of closely associated metabolites (Supplemental Figure 1A). After adjusting for age and gender, two modules differed between MS patients and HCs: a “green” module enriched in γ-glutamyl and oxidative stress metabolites and a “brown” module enriched in benzoate and xanthine metabolites (Figure 3A). We detected some evidence of correlation between metabolites belonging to either module, suggesting that the two distinct modules may be joint predictors of diseases status (Figure 3B; r = 0.27; P = 0.045).

Figure 3

Cross-sectional cohort — results of pathway analyses using weighted correlation network analysis. (A) Box plots of eigen-metabolite levels for green module and brown module that significantly differ between multiple sclerosis (MS) cases and healthy controls (HC) (P = 0.05 for both; based on linear models adjusted for age and sex). Bounds of the box itself represent interquartile range (IQR), while the central line within the box represents the median, and minimum and maximum whiskers represent Q1–Q1.5 × IQR (or the absolute maximum, if smaller) and Q3 + 1.5 × IQR (or the absolute minimum, if larger), respectively. (B) Correlations between metabolites belonging to the green and brown modules. Each pie represents one metabolite within the module, where darker or brighter green or brown hues denote more significant results from the individual tests (P values are listed in Table 2) assessing whether a given metabolite significantly differs between MS patients and HCs. The length of each pie is proportional to the number of correlations a given metabolite has with other metabolites in the green or brown module that exceed 0.45. A connecting line is drawn between the metabolite pies if the metabolite-metabolite correlation is at least 0.45. For example, weighted correlation network analysis (WGCNA) classifies the metabolite glutamate into the green module, and the correlation between maleate and glutamate (classified by WGCNA into the brown module) is ≥0.45. As a result, there is a line connecting the glutamate and maleate pies within the circos plot. We observe some potential intermodular correlation between brown and green modules, where green and brown eigen-metabolites are marginally correlated (Pearson’s r = 0.27; P = 0.045). The contents of these modules are listed in Table 3 and Supplemental Table 1. (C) Relation between metabolite module-membership scores and difference in mean for metabolites between MS cases and HCs. The hue of brown or green color denotes the degree of significance for a test of differences between mean metabolite levels between MS patients and HCs (the darker the hue of green or brown denotes the more significantly different a given metabolite is between MS patients and HCs). Metabolite module-membership scores are derived as the correlation between the overall metabolite module eigen-metabolite score and that of the individual metabolite.

In examining the contents of both modules, we detected significant correlations between module-membership scores (the correlation between module eigen-metabolite level and individual metabolite level is used to estimate the relative importance of a metabolite in a module, a score indicative of the centrality of a metabolite within the network of the module) and the standardized difference between MS patients and HCs. Here, higher metabolite module-membership scores correlated strongly with higher standardized mean differences between MS patients and HCs (Figure 3C; for green module r = 0.67; P < 0.001; for brown module r = 0.61; P = 0.002). Notably, within the green module, levels of γ-glutamyl leucine and γ-glutamyl valine were both correlated with overall eigen-metabolite level and were also highly statistically different (P = 0.002 and P = 0.003, respectively) between MS patients and HCs (Table 3 and Supplemental Table 1; supplemental material available online with this article; https://doi.org/10.1172/jci.insight.95302DS1).

Table 3

Cross-sectional pathway analyses using WGCNA

In complementary analyses, in which we grouped related metabolites based on biological class/function prior to performing any analysis (WGCNA is agnostic), we observed consistent results; metabolites involved in γ-glutamyl amino acid metabolism were significantly altered in MS patients versus HCs (P = 0.04). In addition, lysine metabolites and xenobiotic sulfate chemicals, both of which contribute metabolites to both the green and brown modules in WGCNA, were also marginally different between MS patients and HCs (lysine metabolism: P = 0.07; xenobiotic sulfate chemicals: P = 0.10).

Several metabolites in the green (ornithine, γ-glutamyl valine, γ-glutamyl isoleucine, cysteine glutathione disulfide, 2-aminoadipate) and brown (3-methyl-catechol sulfate, ferulic acid-4-sulfate, N2-furoyl-glycine, caffeine) modules also overlapped with those previously identified in the PLS-DA analysis.

Longitudinal cohort of MS patients and HCs receiving vitamin D supplementation

Pathway analyses using WGCNA and a priori–defined pathway sets. In the separate cohort of MS patients and HCs from the longitudinal vitamin D supplementation study, metabolomic profiling detected 613 compounds, of which 529 (86%) passed quality control procedures and were eligible for analyses. The demographic and disease characteristics of the longitudinal cohort are listed in Table 1.

We again utilized WGCNA to cluster metabolites and identify modules of highly related metabolites. We correlated each of the 14 distinct modules that we extracted with baseline clinical covariates, including age, BMI, and baseline vitamin D status. Baseline 25(OH)D status was significantly correlated with one module eigen-metabolite level (r = 0.42; P = 0.002; Supplemental Figure 2B).

In models assessing longitudinal changes in metabolite modules after vitamin D supplementation, we detected significant differences between MS patients and HCs for two modules: the green module (enriched in γ-glutamyl amino acids and oxidative stress–related metabolites; P for interaction = 0.006; Figure 4A) and the red module (enriched in lysolipids and fatty acids; P for interaction = 0.03; Figure 4B). We did not detect evidence of intermodule correlation between changes in these modules, suggesting that the differential changes in module levels due to vitamin D supplementation are independent (r = 0.03; P = 0.92; Figure 4C).

Figure 4

Vitamin D supplementation cohort — results of pathway analyses using weighted correlation network analysis. Difference in average green metabolite module (A) and red metabolite module (B) changes occurring between multiple sclerosis (MS) cases and healthy controls (HCs). P values are derived from tests of 3-way cross-product between vitamin D status, time, and disease status in generalized estimating equations (GEE) models. (C) Correlations between metabolites belonging to the green and red modules. Each pie represents one metabolite within each of these modules, where darker or brighter green or red hues denote how strongly metabolites changed over time between MS patients and HCs following vitamin D supplementation. The length of each pie is proportional to the number of correlations a given metabolite has with other metabolites in the green or red modules that exceed 0.45. A connecting line is drawn between the metabolite pies if the metabolite-metabolite correlation is at least 0.45. For example, weighted correlation network analysis (WGCNA) classifies the metabolite γ-glutamyl histidine into the green module, and its correlation with γ-glutamyl valine (classified by WGCNA into the green module) is ≥0.45. As a result, there is a line connecting the γ-glutamyl valine and γ-glutamyl histidine pies. In contrast to the modules derived in the cross-sectional cohort (B), we observe no intermodular correlation between green and red modules (Pearson’s r = 0.03; P = 0.71). The contents of these modules are listed in Table 4 and Supplemental Table 2. (D) Relation between metabolite module-membership scores and difference in change in means between MS patients and HCs for the individual metabolite models. The hue of green or red color denotes the degree of significance for a test of differences in the change between MS patients and HCs occurring after vitamin D supplementation for a given metabolite. The darker the red or green hue, the more strongly the metabolite changes differently between MS patients and HCs following vitamin D supplementation. Metabolite module-membership scores are derived as the correlation between the overall metabolite module eigen-metabolite score and that of the individual metabolite.

In analyses of the module contents, there were significant correlations between the eigen-metabolite module score and the effect estimates from individual metabolite models. Higher membership scores corresponded to higher effect estimates (Figure 4D; for the green module: r = 0.43; for the red module: r = 0.61). Within the green module, changes in γ-glutamyl leucine (P = 0.01), γ-glutamyl histidine (P < 0.0001), γ-glutamyl isoleucine (P = 0.002), and γ-glutamyl valine (P = 0.0007) were significantly altered between MS patients and HCs over time. They also appear most strongly related to the module’s eigen-metabolite level, suggesting that these metabolites contribute most to the overall metabolite-module association (Table 4 and Supplemental Table 2). Of the metabolites in the red module with high module-membership scores, significant differential changes in 1-palmitoyl-glycerophosphatidyl inositol (GPI) (16:0), 1-linoleoyl-glycerophosphatidyl choline (18:2), and 1-oleoyl-GPI (18:1) occurred after supplementation. Similarly, in analyses specifying metabolite sets a priori, we observed consistent results; the greatest difference in the rate of change in metabolite levels between MS patients and HCs was for those involved in γ-glutamyl amino acid metabolism (P < 0.001) and lysolipid metabolism (P = 0.07).

Table 4

Vitamin D cohort pathway analyses

Our results show the ability of untargeted metabolomics to identify alterations in the metabolome of MS patients compared with HCs. In particular, we note differences in metabolites related to oxidative stress. Our longitudinal metabolomics study then suggests a differential effect of vitamin D supplementation on the metabolome in MS patients and HCs, with a minimal effect of vitamin D on lowering oxidative stress markers in MS patients. We utilized both a priori and agnostic approaches and noted a consensus among the findings of all statistical approaches employed.

Consistent findings across the analyses implicate altered oxidative stress metabolism in MS patients. Each method identified sets of metabolites related to impaired redox homeostasis and glutathione metabolism (e.g., changes in γ-glutamyl amino acids) (14). These findings are consistent with a large body of literature implicating a role of oxidative stress in MS (15, 16). Elevated oxidative stress markers can be detected in body fluids from patients with MS, and CNS tissue from MS patients displays evidence of oxidative stress–related damage (17, 18). Increased oxidative stress may also play a role in MS progression (16). Therapies targeting oxidative stress by increasing reduced glutathione levels have demonstrated beneficial effects in experimental autoimmune encephalomyelitis (EAE) (19).

Another consistent finding across statistical methods was altered amino acid metabolism. In the PLS-DA model, changes were noted in histidine, phenylalanine, and lysine metabolism. In addition to γ-glutamyl metabolites, WGCNA also identified multiple branched chain amino acids and phenylalanine and lysine metabolites that were altered in MS patients. Metabolomics studies in rodent EAE have demonstrated changes in branched chain amino acid and histidine metabolism (20, 21). Previous studies in MS patients have also noted differences in amino acid metabolism (branched-chain amino acid, alanine, arginine) (22, 23) or in multiple amino acids (tryptophan, alanine, and isoleucine) in addition to energy metabolites (24). Elevations in lactate, certain fatty acids, and N-acetyl species and reductions in glucose and phosphocholine as well as differences in phospholipid metabolism have also been observed in MS (6, 8). We also noted the presence of multiple N-acetyl compounds in the two modules that differed between MS and HC groups and in metabolites identified by PLS-DA (Tables 2 and 3 and Supplemental Table 1). The discrepancies in metabolites identified by these studies could relate to differences in sample handling and processing, analytical methods, and small sample sizes.

Several xenobiotic metabolites contributing to the brown module that significantly differed between MS patients and controls were reduced in the MS group. Levels of benzoate (produced by bacterial metabolism of phenylalanine and polyphenols) and its metabolites were reduced in MS patients. Other xenobiotic compounds also differed in MS patients (e.g., coffee/xanthine metabolites), which could reflect changes in metabolism of the gut microbiota, a recently identified contributor to MS (25, 26). However, due to the cross-sectional design of this study, reverse causation cannot be eliminated.

This study also suggests that metabolomics may be useful in monitoring the effects of therapeutic interventions. We demonstrate that vitamin D supplementation altered two metabolite modules differently in the MS and HC groups. The first module that changed with vitamin D supplementation contained metabolites involved in redox homeostasis: glutathione metabolism and protein oxidation. This is consistent with current knowledge about the molecular effects of vitamin D on related pathways, such as the Nrf2 antioxidant pathway, which leads to activation of the glutathione-S-transferase A2 enzyme, and the antiaging protein Klotho (27, 28). This is also consistent with the relationship between vitamin D deficiency and elevated markers of oxidative stress and the reduction of markers of oxidative stress with vitamin D supplementation in multiple populations (12, 13, 29).

The second module that changed with vitamin D supplementation contained lysolipids and fatty acids. Lysolipids are phospholipids that have had one or more acyl group removed (by hydrolysis). Vitamin D has multiple effects on lipid metabolism through its actions on nuclear hormone receptors, such as the vitamin D receptor and PPARγ (30), including changing the phospholipid content of cells (31). Thus, our results are consistent with the effects of vitamin D on lipid metabolism mediated through these receptors.

Intriguingly, we noted that the effects of vitamin D supplementation on the metabolome were blunted in the MS group compared with the HC group, even after adjusting for differences in the change in serum 25(OH)D levels and other covariates. This suggests a resistance to the metabolic effects of vitamin D in MS patients, which could relate to alterations in vitamin D metabolism or signaling. In trials of vitamin D supplementation in other populations, the effects on oxidative stress and inflammation were influenced by the VDR genotype (32, 33). Thus, the presence of polymorphisms affecting vitamin D metabolism or receptors, or vitamin D–dependent interactions of risk variants and ongoing inflammation, could potentially explain the blunted effects of vitamin D supplementation on the metabolome. An effect of inflammation on the metabolism of vitamin D in MS patients could potentially also explain our findings. However, these hypotheses remain to be confirmed.

The results of our study do have noteworthy limitations. For example, in the cross-sectional study, by design, we could not determine whether the noted differences are a contributing cause to MS or occur as a result of MS. Further longitudinal studies are necessary to properly address this concern. In addition, our study included only people with RRMS, and we cannot comment on potential differences in metabolomics for progressive patients. Another limitation of the study was the lack of data regarding time of blood draw and last meal, which can affect certain metabolites. Limitations of the vitamin D supplementation study include that the cohort focused on women of European descent with low vitamin D levels, limiting the generalizability of our findings. The effects of DMTs on the metabolome are not known; most patients in our study were untreated or were on glatiramer acetate, and, while this is a limitation of the current study, since most MS patients are currently treated with DMTs, it is likely that future studies will also need to contend with this issue. Further research to clarify the influence of genotype on the response of the metabolome to vitamin D supplementation would also be useful.

Using a global untargeted metabolomics approach, we have demonstrated that MS patients appear to have changes in oxidative stress capacity and that vitamin D supplementation reduces markers of oxidative stress and lipid metabolism in HCs but not in MS patients, suggesting an attenuated effect of vitamin D supplementation on the metabolome in MS patients, particularly for metabolites related to oxidative stress. Multiple clinical trials are currently assessing the efficacy of vitamin D supplementation in reducing MS disease activity (34, 35). Understanding the mechanism underlying the reduced effect of vitamin D supplementation on oxidative stress in MS is an important step toward optimizing vitamin D’s therapeutic benefit. Further, since heightened oxidative stress was the most important change in the metabolome in MS patients, interventions that target this finding could result in prognostic improvements. Metabolomics may thus be useful in future trials to monitor the effects of other pharmacological and nonpharmacological therapies on markers of oxidative stress or other relevant pathways.

Participants

For the cross-sectional cohort, patients with RRMS (n = 27) and age-, sex-, and race-matched HCs (n = 27) were recruited from the Johns Hopkins Multiple Sclerosis Center by convenience sampling. Participants provided blood samples (to obtain plasma) and information regarding demographic and disease characteristics. For the vitamin D supplementation study, detailed information on eligibility criteria, recruitment, and follow-up of participants has been previously published (11). Briefly, women of European descent, aged 18 to 60 years old, with BMIs in the range of 18–30 kg/m² and vitamin D levels ≤30 ng/ml, who either had RRMS or were HCs, received 5,000 IU vitamin D₃ daily for 90 days. Individuals with MS were treatment stable and did not experience any attacks over the course of the study of the vitamin D treatment. Figure 1 depicts the study design.

Metabolomics analysis

Samples were stored at –80°C until processed. Sample preparation and mass spectrometry analyses were carried out, as described previously, at Metabolon Inc. (see the Supplemental Methods for details and refs. 36, 37) in separate runs for the cross-sectional and longitudinal cohort. Briefly, recovery standards were added prior to the first step in the extraction process for quality control purposes. To remove protein, to dissociate small molecules bound to protein or trapped in the precipitated protein matrix, and to recover chemically diverse metabolites, proteins were precipitated with methanol under vigorous shaking for 2 minutes (Glen Mills Genogrinder 2000), followed by centrifugation. The resulting extract was divided into 5 fractions: one for analysis by ultra-high performance liquid chromatography–tandem mass spectrometry (UPLC-MS/MS; positive ionization), one for analysis by UPLC-MS/MS (negative ionization), one for the UPLC-MS/MS polar platform (negative ionization), one for analysis by gas chromatography–mass spectrometry, and one sample was reserved for backup (36, 37). Metabolite identification was performed by automated comparison of the ion features in the study samples to a reference library of standard metabolites. Quantification of peaks was performed using area under the curve. Raw values for the area counts for each metabolite were normalized (correcting for variation resulting from instrument interday tuning differences) by the median value for each run day; therefore, medians are 1.0 for each run. We included blinded duplicates (n = 6 pairs) to assess metabolite reliability; median intraclass correlation coefficient for detected metabolites was 0.77 (interquartile range: 0.53–0.91).

The metabolites that were detected by the platform included superpathway (subpathways) amino acids (glycine, serine and threonine, alanine and aspartate, glutamate, histidine, lysine, phenylalanine and tyrosine, tryptophan, leucine, isoleucine and valine, methionine, cysteine and taurine, creatine, urea cycle, arginine and proline, polyamine, guanidine and acetamido, and glutathione metabolism), peptides (γ-glutamyl amino acid, dipeptide derivative, dipeptide, polypeptide, fibrinogen cleavage peptide, acetylated peptide), carbohydrates (glycolysis, gluconeogenesis and pyruvate, pentose, glycogen, disaccharides and polysaccharides, fructose, mannose and galactose, amino-sugar, advanced glycation end products), energy metabolites (TCA cycle and oxidative phosphorylation), lipids (short-chain fatty acid, medium-chain fatty acid, long-chain fatty acid, PUFA, branched fatty acid, dicarboxylate fatty acid, amino fatty acid, fatty acid metabolism, acyl glycine, acyl carnitine, carnitine, ketone bodies, monohydroxy fatty acid, dihydroxy fatty acid, oxidized Fa, endocannabinoid, eicosanoid, inositol, phospholipid, phosphatidylserine, lysolipid, plasmalogen, lysoplasmalogen, glycerolipid, monoacylglycerol, diacylglycerol, sphingolipid, mevalonate, sterol, steroid, primary bile acid, secondary bile acid, acyl choline, acyl glutamine), nucleotides (purine, pyrimidine), cofactors and vitamins (tocopherol, riboflavin, pantothenate and CoA, ascorbate and alderate, nicotinamide and nicotinate, biotin, hemoglobin and porphyrin, vitamin A, vitamin B6 metabolism), and xenobiotics (benzoate, tobacco, xanthine, food component, bacterial/fungal, drug, chemical).

Vitamin D level measurement

Serum 25(OH)D levels for the vitamin D supplementation cohort were measured in a single batch at Heartland Assays using liquid chromatography–mass spectrometry.

Statistical analysis

Cross-sectional study. Metabolite concentrations obtained from the previous steps underwent preprocessing, including removal of metabolites with greater than 30% missing values. We imputed missing values using the k-nearest neighbors algorithm (consistent results were observed when missing values were imputed as the minimum of observed metabolites). We then applied log transformations and standardized the resultant log-transformed values. We applied principal components analysis to identify outliers, followed by PLS-DA to assess the ability of the metabolites to differentiate MS patients from HCs. To avoid overfitting/overinterpreting results of PLS-DA, we performed cross-validation of model-validity measures by permuting the class labels (1,000 permutations) and generating P values for the R²Y and Q²Y statistics based on this distribution; this allowed us to identify outliers based on orthogonal distance from the model plane. Following assessment of model validity, VIP values were obtained for all metabolites to identify the metabolites that drive separation of the two groups.

Because our study sample size was relatively small in relation to the number of metabolites quantified, we focused our analytic approach on assessing metabolic pathway/set-wide differences between MS patients and HCs. We performed two sets of complementary analyses: an agnostic approach to discover pathways and another incorporating known information on metabolic pathways. For the agnostic approach, we applied WGCNA to identify related metabolite modules (38). WGCNA is a systems biology method that was originally developed to characterize correlation patterns among genes from microarray studies. It has been applied in other settings, such as metabolomics, cancer, and analysis of brain imaging data (39–42). As an analytic technique, its primary goal is to find clusters (modules) of interconnected nodes (e.g., metabolites in this case), and it allows summary measures of the clusters to be used in subsequent analysis. This helps to alleviate some of the multiple testing problems introduced by analyzing hundreds of nodes because module scores are used in the analysis rather than an individual marker. WGCNA can also identify highly connected hub nodes (metabolites), which are centrally located within the module, and prior studies of gene expression have shown that such hub nodes are more likely to be biologically relevant markers. Briefly, WGCNA identifies modules by constructing a correlation network of metabolites and derives eigen-metabolite scores (corresponding to the first principal component) from the identified modules.

These scores were used to test for a differential association between each module and MS patients versus HCs (adjusting for the matching factors age and sex using unmatched generalized linear models). We also identified potential hub nodes or metabolites within each module with high intramodular importance (as shown through a strong correlation between eigen-metabolite score for a given module with the individual metabolite level). Metabolite modules were derived among the HCs to ensure that the distribution of correlations between metabolites was not distorted by disease status.

For the a priori– defined pathway analyses, we classified metabolites into groups (≥5 metabolites) based on biologic function (e.g., glutathione metabolism; xanthine metabolism) and applied a resampling permutation-based algorithm to assign statistical significance while preserving metabolite-metabolite correlation. This approach is similar to permutation based approaches adopted for pathway analysis of gene expression data. Briefly, we fit individual models for each metabolite and extracted and ranked the P value, which was used to calculate the average rank of P values for a given pathway set. We then permuted phenotype labels 1,000 times and repeated the above procedure to calculate an average rank of the P values for a pathway set for each permutation. The P value for a given pathway set is the probability the observed average rank is less than the expected average rank (derived from the 1,000 permutations).

Vitamin D supplementation study. We applied similar quality control procedures and preprocessing, as in the cross-sectional cohort. Since the study was performed at two sites (Johns Hopkins and UCSF), we removed any batch-type effects of site using the Combat algorithm (43). We performed analyses based on metabolite pathway-based sets using WGCNA and a priori–defined metabolite pathway sets.

After deriving metabolic modules using WGCNA, we used the eigen-metabolite from each module as the dependent variable in regression models, including clinical traits, time, and vitamin D status. We tested whether changes in eigen-metabolite levels associated with changes in vitamin D status differed between MS patients and HCs using generalized estimating equations (GEE) and a 3-way interaction term (the cross-product between case status [MS or HC], time [baseline or follow-up], and vitamin D status) and including all lower order interaction terms. For eigen-metabolite modules that differed between MS cases and HCs, we also assessed both intramodule and intermodule relatedness using similar procedures as in the cross-sectional cohort. For the a priori–defined pathway analyses, we fit similar GEE models (as in the eigen-metabolite analyses), but using the individual metabolite level as the dependent variable, and extracted the P value for the cross-product 3-way interaction term (MS or HC, time, vitamin D status)and subsequently ranked these P values. We calculated the average rank for each of the a priori–defined pathway sets and assigned statistical significance using resampling-based methods as the cross-sectional analyses. All models were additionally adjusted for age and BMI.

P values of less than 0.05 were considered statistically significant.

Study approval

Both study protocols were approved by the institutional review boards of Johns Hopkins University, Baltimore (for the cross-sectional and longitudinal cohort), and University of California, San Francisco (for the longitudinal cohort). All participants provided informed consent prior to any study procedures being performed.