Changes in gut microbiota have been associated with several diseases. Here the international Multiple Sclerosis Microbiome Study (iMSMS) studied the gut microbiome of 576 MS patients (36% untreated), and genetically unrelated household healthy controls (1,152 total subjects). We observed a significantly increased proportion of Akkermansia muciniphila, Ruthenibacterium lactatiformans, Hungatella hathewayi, and Eisenbergiella tayi and decreased Faecalibacterium prausnitzii and Blautia species. The phytate degradation pathway was over-represented in untreated MS, while pyruvate-producing carbohydrate metabolism pathways were significantly reduced. Microbiome composition, function and derived metabolites also differed in response to disease-modifying treatments. The therapeutic activity of interferon-β may in part be associated to upregulation of short-chain fatty acid transporters. Distinct microbial networks were observed in untreated MS and healthy controls. These results strongly support specific gut microbiome associations with MS risk, course and progression, and functional changes in response to treatment.

Recruitment and inclusion criteria

A total of 576 MS patients and their household healthy controls were included in this study. The first 128 MS-control pairs were recruited as Cohort 1 (Zhou et al., 2020) and the subsequent 448 pairs were recruited as Cohort 2. Participants were recruited through MS clinics at UCSF (San Francisco, CA), Brigham and Women’s Hospital (Boston, MA), Mount Sinai (New York, NY), the Anne Rowling Clinic (Edinburgh, UK), University of Pittsburgh (Pittsburgh, PA), Biodonostia Health Research Institute (San Sebastián, Spain) and FLENI (Buenos Aires, Argentina). Each collaborating site obtained human subject research approval through their respective ethics review committees, following a master protocol established at UCSF (protocol no. 15-17061). All participants provided written informed consent and signed a HIPAA Authorization allowing for the use of their medical record for research purposes. Each collaborating site obtained human subject research approval through their respective ethics review committees, following a master protocol established at UCSF (protocol no. 15-17061). All participants provided written informed consent and signed a HIPAA Authorization allowing for the use of their medical record for research purposes.

Inclusion criteria required that participants carry a diagnosis of MS;(McDonald et al., 2001) be of White (Hispanic or non-Hispanic) ethnicity (i.e. to match characteristic genetic risk profile of MS(Baranzini and Oksenberg, 2017)); and be enrolled with a genetically unrelated household control with cohabitation for at least six months. Exclusion criteria for MS and control subjects included the presence of other autoimmune disorders, gastrointestinal infections, and other neurological disorders. Participants were excluded if they received oral antibiotics within the past three months, corticosteroids within the past 30 days, or were on a disease-modifying therapy for less than three months.

Specimen collection

Participants were provided with a stool sample collection kit, and instructed to obtain two consecutive stool samples in the privacy of their own homes. Each stool sample time point included three collection vials—a Q-tip (Q, dry), a snap-frozen vial (S, wet), and a vial filled with lysogeny broth (LB) and glycerol. Participants were instructed to freeze the samples for at least 12 hours, and ship them frozen. Samples were returned to each site via overnight shipping. Blood samples were collected at the initial visit only and stored at −80°C upon further processing. Blood samples were collected at the initial visit only and stored at −80°C upon further processing.

Stool sample preparation for sequencing

For the first cohort, Q-tip samples (i.e. dry) and snap frozen (i.e. wet) samples were processed using the QIAamp PowerFecal DNA Kit (ref 12830-50).  After lysis solution was added to bead beating tubes, dry samples were transferred by grinding the Q-tips into the bottom while snap frozen samples were chipped to an appropriate size for the kit.  Sample processing was done on a QIAcube platform according to the protocols generated by the manufacturer (QIAGEN).  DNA sample quantity and purity were measured by NanoDrop spectrophotometry (Thermo Scientific^TM).  The second cohort samples were processed using the MagAttract PowerSoil DNA EP Kit (ref 27100-4-EP).  After lysis solution was added to the bead beating plate, samples were added to each well in the same manner used previously for bead beating tubes.  Physical lysis was executed using a mixer mill and subsequent steps were automated using the EpMotion platform. Sample quality and quantity were assessed with the same method used for the first cohort.

16S rRNA sequencing

The V4 region of the bacteria 16S ribosomal RNA gene was amplified on an Illumina MiSeq platform using the Earth Microbiome Project protocol (Caporaso et al., 2012). Amplicon reads from two cohort samples were analyzed using QIITA(Gonzalez et al., 2018; Hillmann et al., 2018) to combine the forward and reverse reads, trim short reads of less than 150bp and assign filtered reads to amplicon sequencing reads (ASVs) using default Deblur parameters against Greengenes (version 13.8 at 99% identity) as described in QIIME2 documents.(Caporaso et al., 2010) As the impact of sample collection method on microbial composition is negligible (Zhou et al., 2020), sequencing counts of samples from each participant were summed. ASVs were filtered to retain only the ones covering at least 10 total reads and present in at least 5% of samples for downstream analyses (Table S4).

Microbial diversity

The ASVs characterized by 16s rRNA sequencing were rarefied to 10,000 sequences per participant sample for microbial diversity analysis. a-diversity was measured by Shannon(Shannon, 1997) and Chao1(Chao, 1984) indexes (Table S4). Both weighted and unweighted UniFrac(Lozupone and Knight, 2005) distances were computed between all samples (Table S4), and principal coordinates analysis (PCoA) was applied to visualize the b-diversity. All these analyses were performed with QIIME2. Bray-Curtis(Bray JR, 1957) dissimilarities were calculated to compare gut microbiome among individuals in terms of geographic distance. Statistical significance was determined by ANOVA. The PERMANOVA test (McArdle and Anderson, 2001) was used to assess the effect of host metadata categories (confounders): demography, lifestyle, diseases, medication, and physiology, on the variation of microbiome abundance (Table S4). The test was performed by using the “adonis” function implemented in R package vegan(Zapala and Schork, 2006) and tested on weighted UniFrac distances of paired MS and HHC samples with reported host factors. The variance of microbial abundance between MS and control or between treated/untreated MS and controls was tested by specifying “strata” as household to control the within-house comparison. The empirical p-value was obtained by running 999 permutations. When appropriate, statistical p-values were adjusted by false discovery rate (FDR).

Shallow whole metagenome shotgun sequencing (WMGS) and data processing

1 ng of input DNA was used in a 1:10 miniaturized Kapa HyperPlus protocol. For samples with less than 1 ng DNA, a maximum volume of 3.5 μl input was used. Library concentration was determined with triplicate readings of the Kapa Illumina Library Quantification Kit (cohort 1, ref 07962428001) or Pico Green Quantification Kit (cohort2, ref P11496); 20 fmol of sample libraries were pooled and size-selected for fragments between 300 and 800 bp on the Sage Science PippinHT to exclude primer dimers. The pooled library was sequenced as a paired-end 150-cycle run on an Illumina HiSeq2500 v2 (cohort1) or NovaSeq 6000 (cohort2) at the UCSD IGM Genomics Center with sequencing depth of 0.5 million reads per sample.

Demultiplexed shallow shotgun metagenomic sequences were processed using Atropos (v1.1.24 ) (Didion et al., 2017) to remove adapters (forward “GATCGGAAGAGCACACGTCTGAACTCCAGTCAC” , reverse “GATCGGAAGAGCGTCGTGTAGGGAAAGGAGTGT”) and filter reads with lower quality score than 15 and length less than 80 base pairs. For taxonomic assignment reads were aligned to the Web of Life(Zhu et al., 2019) of 10,575 bacterial and archaeal genomes using SHOGUN (Hillmann et al., 2018) in the Bowtie2 alignment mode. Species-level functional profiling was performed using HUMAnN2 default. (Franzosa et al., 2018) Sequencing counts of samples from each participant were summed (Table S5).

Microbial co-abundance network

Co-abundance network inference was performed on WMGS species using SparCC(Friedman and Alm, 2012) method (in R using SpiecEasi package(Kurtz et al., 2015)), which is a tool to infer linear relationships with high precision for high diversity compositional data. SparCC correlations were adjusted for age, sex, and BMI using cor2por function from R package “corpcor”. Significant co-abundance was controlled at FDR 0.05 level using 100 × permutation (Table S6). In each permutation, the abundance of each microbial factor was randomly shuffled across samples. To keep the co-abundances with high correlations in a dense microbial network, we filtered co-abundances with a lower absolute correlation than 0.4 and subnetworks with only two species.

Metabolite measurement in stool and serum samples

Blood samples were centrifuged at 2200g for 20 minutes. The serum layers were aspirated and moved into 2mL cryovials. The serum samples were stored at −80°C before metabolomics profiling. Fecal (150g/sample) and serum (150ul/sample) samples were shipped on dry ice to Metabolon Inc. (Durham, North Carolina) and maintained at −80°C until processed following their published protocols.(Evans et al., 2009; Long et al., 2017; McMurdie et al., 2022) Global metabolite and 8 targeted short-chain fatty acid profiling were analyzed in a subset of samples (n = 490), including untreated RRMS patients (N=79), and in those treated with dimethyl fumarate (n=47), fingolimod (n=39), glatiramer acetate (n=31), and interferon-β (n=49), as well as their corresponding household healthy controls (Table S7).

Global metabolomic profiling

Samples were prepared using the automated MicroLab STAR® system from Hamilton Company.  Several recovery standards were added prior to the first step in the extraction process for QC purposes. To remove protein, dissociate small molecules bound to protein or trapped in the precipitated protein matrix, and recover chemically diverse metabolites, proteins were precipitated with methanol under vigorous shaking for 2 min (Glen Mills GenoGrinder 2000) followed by centrifugation.  The resulting extract was divided into five fractions: two for analysis by two separate reverse phases (RP)/UPLC-MS/MS methods with positive ion mode electrospray ionization (ESI), one for analysis by RP/UPLC-MS/MS with negative ion mode ESI, one for analysis by HILIC/UPLC-MS/MS with negative ion mode ESI, and one sample was reserved for backup.  Samples were placed briefly on a TurboVap® (Zymark) to remove the organic solvent.  The sample extracts were stored overnight under nitrogen before preparation for analysis. 

Several types of controls were analyzed in concert with the experimental samples: a pooled matrix sample generated by taking a small volume of each experimental sample (or alternatively, use of a pool of well-characterized human plasma) served as a technical replicate throughout the data set; extracted water samples served as process blanks; and a cocktail of QC standards that were carefully chosen not to interfere with the measurement of endogenous compounds was spiked into every analyzed sample, allowed instrument performance monitoring and aided chromatographic alignment.  Instrument variability was determined by calculating the median relative standard deviation (RSD) for the standards that were added to each sample prior to injection into the mass spectrometers.  Overall process variability was determined by calculating the median RSD for all endogenous metabolites (i.e., non-instrument standards) present in 100% of the pooled matrix samples.  Experimental samples were randomized across the platform run with QC samples spaced evenly among the injections.

All methods utilized a Waters ACQUITY ultra-performance liquid chromatography (UPLC) and a Thermo Scientific Q-Exactive high resolution/accurate mass spectrometer interfaced with a heated electrospray ionization (HESI-II) source and Orbitrap mass analyzer operated at 35,000 mass resolution.  The sample extract was dried and then reconstituted in solvents compatible to each of the four methods. Each reconstitution solvent contained a series of standards at fixed concentrations to ensure injection and chromatographic consistency.  One aliquot was analyzed using acidic positive ion conditions, chromatographically optimized for more hydrophilic compounds. In this method, the extract was gradient eluted from a C18 column (Waters UPLC BEH C18-2.1x100 mm, 1.7 µm) using water and methanol, containing 0.05% perfluoropentanoic acid (PFPA) and 0.1% formic acid (FA). Another aliquot was also analyzed using acidic positive ion conditions, however, it was chromatographically optimized for more hydrophobic compounds.  In this method, the extract was gradient eluted from the same aforementioned C18 column using methanol, acetonitrile, water, 0.05% PFPA and 0.01% FA and was operated at an overall higher organic content. Another aliquot was analyzed using basic negative ion optimized conditions using a separate dedicated C18 column. The basic extracts were gradient eluted from the column using methanol and water, however with 6.5mM Ammonium Bicarbonate at pH 8. The fourth aliquot was analyzed via negative ionization following elution from a HILIC column (Waters UPLC BEH Amide 2.1x150 mm, 1.7 µm) using a gradient consisting of water and acetonitrile with 10mM Ammonium Formate, pH 10.8. The MS analysis alternated between MS and data-dependent MSn scans using dynamic exclusion.  The scan range varied slightly between methods but covered 70-1000 m/z. Raw data files are archived and extracted as described below.

Raw data were extracted, peak-identified and QC processed using Metabolon’s hardware and software. These systems are built on a web-service platform utilizing Microsoft’s .NET technologies, which run on high-performance application servers and fiber-channel storage arrays in clusters to provide active failover and load-balancing. Compounds were identified by comparison to library entries of purified standards or recurrent unknown entities. Metabolon maintains a library based on authenticated standards that contains the retention time/index (RI), mass to charge ratio (m/z), and chromatographic data (including MS/MS spectral data) on all molecules present in the library. Furthermore, biochemical identifications are based on three criteria: retention index within a narrow RI window of the proposed identification, accurate mass match to the library +/- 10 ppm, and the MS/MS forward and reverse scores between the experimental data and authentic standards. The MS/MS scores are based on a comparison of the ions present in the experimental spectrum to the ions present in the library spectrum. While there may be similarities between these molecules based on one of these factors, the use of all three data points can be utilized to distinguish and differentiate biochemicals. More than 3300 commercially available purified standard compounds have been acquired and registered into LIMS for analysis on all platforms for determination of their analytical characteristics.  Additional mass spectral entries have been created for structurally unnamed biochemicals, which have been identified by virtue of their recurrent nature (both chromatographic and mass spectral). These compounds have the potential to be identified by future acquisition of a matching purified standard or by classical structural analysis.

A variety of curation procedures were carried out to ensure that a high-quality data set was made available for statistical analysis and data interpretation. The QC and curation processes were designed to ensure accurate and consistent identification of true chemical entities and to remove those representing system artifacts, misassignments, and background noise.  Metabolon data analysts use proprietary visualization and interpretation software to confirm the consistency of peak identification among the various samples. Library matches for each compound were checked for each sample and corrected if necessary.

Targeted short-chain fatty acid profiling

Human feces and human serum samples are analyzed for eight short chain fatty acids: acetic acid (C2), propionic acid (C3), isobutyric acid (C4), butyric acid (C4), 2-methyl- butyric acid (C5), isovaleric acid (C5), valeric acid (C5), and caproic acid (hexanoic acid, C6), with the addition of lactic acid by request, by LC-MS/MS (Metabolon Method TAM135: “LC-MS/MS Method for the Quantitation of Short Chain Fatty Acid (C2 to C6) in Human Feces” and TAM148: “LC-MS/MS Method for the Quantitation of Short Chain Fatty Acid (C2 to C6) in Human Plasma and Serum”). Human feces and human serum samples are spiked with stable labelled internal standards and are homogenized and subjected to protein precipitation with an organic solvent. After centrifugation, an aliquot of the supernatant is derivatized. The reaction mixture is diluted, and an aliquot is injected onto an Agilent 1290/AB Sciex QTrap 5500 LC MS/MS system equipped with a C18 reversed-phase UHPLC column. The mass spectrometer is operated in negative mode using electrospray ionization (ESI). The peak area of the individual analyte product ions is measured against the peak area of the product ions of the corresponding internal standards. Quantitation is performed using a weighted linear least squares regression analysis generated from fortified calibration standards prepared immediately prior to each run. LC-MS/MS raw data are collected using AB SCIEX software Analyst 1.6.2 and processed with SCIEX OS-MQ software v1.7.

Diet data

A validated Block 2005 food frequency questionnaire (FFQ)(Block, 2005) was set up through an external vendor (NutritionQuest). The intake of foods and nutrients was measured by NutritionQuest in a standardized fashion for all participants based on their responses to the FFQ. 37 nutrient items were summarized and grouped as antioxidants, average intake, B-vitamins, food group servings, and minerals (Table S2). The Healthy eating index (HEI2015 with 10 components) was also calculated for all qualifying participants through NutritionQuest (Table S3).