Gut microbiome of multiple sclerosis patients and paired household healthy controls reveal associations with disease risk and course

of


Introduction
Multiple sclerosis (MS) is an autoimmune disease of the central nervous system (CNS) characterized by demyelination, axonal damage, and progressive neurologic disability. The etiology and pathogenesis of MS is complex and remain elusive, although both genetic and environmental factors are involved. Gut microbiota, an important modulator of the immune response (Geva-Zatorsky et al., 2017) and brain function, has emerged as a likely environmental contributor to MS. (Esmaeil Amini et al., 2020;Kadowaki and Quintana, 2020;Probstel and Baranzini, 2018) Alterations in commensal gut microbiota (referred to as dysbiosis) have been linked to many inflammatory conditions. (Honda and Littman, 2016) Numerous studies including ours have shown both depletion and enrichment of certain bacteria in MS patients compared with healthy controls, (Berer et al., 2017;Cekanaviciute et al., 2017;Chen et al., 2016a;Cox et al., 2021;Jangi et al., 2016) suggesting certain taxa might be associated to either disease pathogenesis or progression. For instance, Akkermansia muciniphila was found to be increased in MS samples and exacerbated intestinal inflammation in mouse models. Jangi et al., 2016) However, a recent report showed a protective effect of MS-derived Akkermansia in experimental allergic encephalomyelitis (EAE), (Cox et al., 2021) suggesting a complex scenario in which context, rather than individual taxa, might be the driving factor of bacteria-human interactions. It remains uncertain as to whether the disease state occurs in response to microbial alterations, or alternatively that the disease drives these gut microbiome alterations. mouse and human studies indicate that microbiota can potentially affect the onset and progression of diseases mediated by different immune effector cells and soluble metabolic, immune, and neuroendocrine factors modulated by gut microbes. (Camara-Lemarroy et al., 2018;Probstel and Baranzini, 2018) Although microbial changes in MS have been detected across different studies, most of the alterations were reported in relapsing-remitting MS (RRMS), whereas few studies investigated the microbiome in progressive MS. Furthermore, it is difficult to identify a common pattern since results are rarely concordant. Kozhieva et al. found an increase of Gemmiger sp. and Ruminococcaceae species in a Russian cohort of 15 primary progressive (PPMS) patients compared to healthy subjects, (Kozhieva et al., 2019) while a Belgian cohort of 26 PPMS patients showed a decrease of Gemmiger and Butyricicoccus. (Reynders et al., 2020) A recent study identified elevated Enterobacteriaceae, Bifidobacterium animalis, Clostridium g24 FCEY, Dorea massiliensis, Longicatena, and Ruminococcaceae FJ366134 in 44 progressive MS patients compared to both RRMS and healthy subjects. (Cox et al., 2021) The gut microbiota can also be altered by drugs with either beneficial or undesirable effects. Many common drugs have antimicrobial effects, or exert a large impact on the composition of gut microbiome. (Maier et al., 2018) Studies in MS and other diseases have suggested that therapy-induced microbial change may contribute to efficacy by favoring microbes with anti-inflammatory properties, or able to metabolize the drug more effectively. (Alexander et al., 2017;Sand et al., 2019) For instance, glatiramer acetate or dimethyl fumarate treatments were associated with an increase in the abundance of common bacterial groups with anti-inflammatory effects, like Bacteroides, in people with RRMS. (Sand et al., 2019) Several RRMS associated microbial changes, (e.g. decrease of Prevotella copri, Faecalibacterium prausnitzii and Roseburia intestinalis), were modulated by disease modifying therapies (DMTs), (Castillo-Alvarez et al., 2018;Cox et al., 2021;Jangi et al., 2016) suggesting that therapeutic efficacy may be due to the effects of DMTs on gut microbiota.
However, current microbiome studies in MS are limited by the relatively small size of the cohort analyzed, and inadequate handling of multiple confounding factors, such as genetic heterogeneity of participants, geographic location, disease subtype, treatment and diet. Also, many studies rely on 16S ribosomal RNA sequencing, which offers low resolution to identify MS associated species. To overcome these challenges, the international MS Microbiome Study (iMSMS) is systematically recruiting MS patients and household healthy controls in the US, Europe and South America. The advantages of the household-controlled experimental design, sequencing method and handling of confounding factors (e.g. geographic location and diet) on gut microbiome were recently reported in a pilot cohort of 128 patient:control pairs.  In this study, we present a large microbiome study of MS and healthy controls (n= 576 pairs), and investigate relationships with MS susceptibility, progression, and treatment.

Recruitment.
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  and the subsequent 448 pairs were recruited as the Cohort 2. Recruitment details have been described elsewhere.  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. 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 cohabitatation 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.

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 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,  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.

Microbial diversity.
The ASVs characterized by 16s rRNA sequencing were rarefied to 10,000 sequences per participant sample for microbial diversity analysis. -diversity was measured by Shannon (Shannon, 1997) and Chao1 (Chao, 1984) indexes. Both weighted and unweighted UniFrac (Lozupone and Knight, 2005) distances were computed between all samples, and principal coordinates analysis (PCoA) was applied to visualize the -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. 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 were 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) or Pico Green Quantification Kit (cohort2); 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 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  of 10,575 bacterial and archaeal genomes using SHOGUN (Hillmann et al., 2018) in the Bowtie2 alignment mode. Sequencing counts of samples from each participant were summed.

Metabolite measurement.
Fecal and serum 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 and analyzed in a sub-set of samples.

Statistical analysis.
Microbial composition was normalized as relative abundance and further transformed with a variance-stabilizing arcsine square-root transformation. (Lloyd-Price et al., 2019;Morgan et al., 2012;Sokal, 1982) Global metabolite intensity and SCFA concentration were normalized by log transformation. Mixed linear regression model was applied on transformed data to identify differential features (species, pathways and metabolites) by adjusting random effects of house and recruitment site, and fixed effects of age, sex and BMI. The linear regression was performed using lmer function from R package "lme4". In order to reduce the effect of zero-inflation in microbiome data, a variance filtering step was applied to remove species features with very low variance (<1E-5). The contribution of individual species in a specific pathway was visualized in a bar plot using HUMAnN2 "humann2_barplot" function. The organism-pathway-reaction-compound network was built by Scalable Precision Medicine Oriented Knowledge Engine (SPOKE), a large graph with multiple types of nodes and relationships integrated from more than 30 publicly available databases covering human and bacterial molecular interactions. (Himmelstein et al., 2017;Nelson et al., 2021) Altered metabolites were linked to gut microbes through reactions (MetaCyc and KEGG) mediated by microbial gene families screened in our WGMS data. Functional KEGG enrichment analysis of metabolites was performed using MetaboAnalyst 5.0. (Pang et al., 2021) To identify species associated with disease severity, the updated global Multiple Sclerosis Severity Score (uGMSSS) was calculated by combining the Expanded Disability Status Scale (EDSS) and disease duration using global_msss function from R package "ms.sev". We focused on the species with prevalence in more than 50% samples, spearman correlations were calculated and tested adjusting for age and BMI using pcor.test function from R package "ppcor".

Microbial co-abundance network.
Co-abundance network inference was performed 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. 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. To test whether the microbial co-abundance relationships showed case or control specificity, i.e. whether the effect size of co-abundance in MS group was very different from that in healthy control, we applied the IQR (interquartile ranges) based the outlier detection method as adapted in paper.  The effect size for co-abundance was measured by the SparCC correlation coefficient in our analysis. The effect sizes were ranked from low to high and extracted corresponding 25%, 50% and 75% quartile values (Q1, Q2 and Q3, respectively). IQR was then calculated as Q3-Q1. The specific co-abundance was defined in each corresponding MS or healthy group if the effect size fell outside of Q1 − 2 × IQR (smallest) or Q3 + 2 × IQR (largest). Cohort specific species were linked to their MetaCyc pathways, as many species share pathways, we focused on those that are unique to the cohort specific species.

Diet analysis.
A validated Block 2005 food frequency questionnaire (FFQ) (Block, 2005) was set up through an external vendor (NutritionQuest). The intake of foods and nutrients were 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. Dietary dissimilarity was measured using Jaccard distance of the nutrient intake. The effect of confounders on the variation of diet and the effect of dietary items (covariates) on the variation of gut microbiome were accessed by PERMANOVA (Permutational multivariate analysis of variance). (McArdle and Anderson, 2001) The test was performed by using the "adonis" function implemented in R package vegan. (Zapala and Schork, 2006) The empirical P-value was obtained by running 999 permutations.

Availability of data and materials.
The dataset generated and analyzed during the current study are available in the EMBL-ENA repository (https://www.ebi. ac.uk) with accession number ERP115476 (Supplementary Data S1). No custom software was written for this manuscript. R code is available at GitHub (https://github.com/BaranziniLab/iMSMS_study).

Results
A total 576 pairs of MS patients and genetically unrelated household healthy controls (HHCs) were recruited between September 2015 and January 2019 from seven sites (recruiting centers) located in San Francisco, Boston, New York, Pittsburgh, Buenos Aires, Edinburgh and San Sebastián ( Figure 1, Table 1, Supplementary Data S1). The first 128 pairs were recruited before October 2016 (Cohort1 ) and the subsequent 448 pairs were recruited before January 2019 (Cohort 2). The majority of participants (97.3%) are Caucasian and/or Hispanic. As expected, an uneven proportion of female MS patients was observed (69.4%), but there was no significant difference in age or body mass index (BMI) between the MS and control groups. The mean Expanded Disability Status Scale (EDSS) score of MS patients was 2.6 (IQR: 1-4) with a mean of disease duration 14.2 years (IQR: 6-21). The median disease duration corrected Multiple Sclerosis Severity Score (MSSS) was 3.37 (IQR: 0.86-5.57). Among the 576 MS patients, 209 (36%) were untreated and 367 (63%) were treated with a disease modifying therapy (DMT). Treatments included oral agents Fingolimod (n=71), and dimethyl fumarate (DMF, n=86);, injectables glatiramer acetate (GA, n=68) and interferon (IFN, n=87); and infusion agents anti-CD20 monoclonal antibodies (n=28), and natalizumab (n=27). Of the total of 576, 437 (76%) patients had RRMS, 68 (12%) secondary progressive MS (SPMS) and 71 (12%) primary progressive (PPMS). Given the heterogeneity in the assessment of patients with SPMS and PPMS, they were combined into a single category, progressive MS (PMS, n=139, 24%) for subsequent analyses. As a natural consequence of the disease process, patients with progressive disease were relatively older, and had higher MSSS than those with RRMS. In addition, 74% RRMS patients were on a DMT, while only 30% of PMS were treated, reflecting the higher effectiveness of available therapies for relapsing disease and the lack of therapies for progressive MS.
All participants completed a clinical survey to report the disease status and treatment, and a high proportion of participants (94%, n = 1,086) completed the subject survey to report the demographics, medication, lifestyle and physiology factors ( Figure 1A, and Supplementary Data S1). The majority of participants (90%, n= 1,034) also completed the online food frequency questionnaire (FFQ). A summary of dietary questionnaires and the dietary intakes, including average intake, food group servings, antioxidants, minerals and vitamins is provided in Supplementary Data S2. The Healthy eating index (HEI2015 with 10 components) was also calculated for all qualifying participants (Supplementary Data S3).

Altered Gut microbial composition in MS
We first used 16s rRNA data to study the global microbial composition (α-diversity and βdiversity). 16S rRNA sequencing has been more commonly used in microbiome studies to date, thus several well-established databases (e.g. Greengenes(DeSantis et al., 2006)) are available. The 576 pairs were processed and sequenced in two cohorts (128 pairs in Cohort 1 and 448 pairs in Cohort 2) (Methods, Supplementary Data S4). The microbial composition and diversity were highly correlated in duplicate samples sequenced in the two cohorts ( Figure S1A-C), and also in duplicated samples processed by different DNA isolation methods ( Figure S1D-F), thus allowing us to merge all sequencing samples for a joint analysis. After removing samples with low coverage (<10,000 reads), 500 pairs of MS and household control participant samples were used for diversity analysis (Supplementary Data S5).
No significant difference in α-diversity was observed between MS and HHC groups as measured by Shannon ( Figure 1B) and Chao1 indexes ( Figure S2A). We also found no significant difference in αdiversity across patient:HHC pairs of RRMS, PMS ( Figure 1B), untreated MS, and treated MS ( Figure S2B). Intriguingly, β-diversity-based sample clustering revealed a significant difference in MS regardless of treatment, and also differed in untreated or treated MS group status compared to their HHC ( Figure 1C, PERMANOVA, p < 0.05). No significant difference was observed between untreated and treated MS patients ( Figure 1C, PERMANOVA, p > 0.05). Different microbial communities were also observed across patient:HHC pairs of RRMS and PMS patients, and when comparing RRMS versus PMS patients ( Figure 1D). We next tested how much of the variance in microbial diversity (weighted uniFrac distances) was explained by the host confounders, including demographics, lifestyle, disease, medication, and physiology factors (Supplementary Data S5). (Vujkovic-Cvijin et al., 2020; Not surprisingly, the recruitment site showed a significant and dominant effect on the microbial composition, as we and other studies have reported ( Figure 1E). (Gaulke and Sharpton, 2018; By checking the gut microbial -diversity in individuals from each recruitment site, we observed lower microbial diversity in both healthy and MS participants from New York whereas a higher diversity in participants from San Francisco and San Sebastián ( Figure S2C). We hypothesize that these differences in microbial diversity might be associated with varied dietary habits (see dietary analysis). Microbial differences associated with geography were also shown by the PCoA of the microbiome β-diversity ( Figure S2D, PERMANOVA, P < 0.05). The second and third largest component was explained by disease status (RRMS/HHC, PMS/HHC) and treatment status (treated MS/HHC, untreated MS/HHC), implying an altered gut microbiome in MS patients versus HHC as well as an effect of treatment on changing microbial structure.
Age, sex and BMI also showed significant effects on microbial compositions ( Figure 1E). Our household design effectively reduces age-associated variation as the great majority of household participants are spouses with comparable age (Table 1). Smoking and education also exerted significant effects, but these were highly different across recruitment sites (e.g. participants from San Francisco, Boston and New York are less likely to smoke and reported higher education) ( Figure S2E-F). Because we detected no significant difference in smoking (Fisher's exact test P > 0.05) and comparable education levels in MS patients and HHC, we hypothesized that these effects would be adjusted by controlling the site effect. A smaller effect was identified by medication use and MS comorbidities ( Figure 1E) as MS patients tend to use more medications and have depression or anxiety ( Figure S2G, Fisher's exact test P < 0.001). No significant microbial divergence was related to factors such as household pets, birth method, or asthma in our study (Supplementary Data S5).

Disease associated microbial changes adjusted for confounders
Whole-metagenome shotgun sequencing (WMGS) provides higher resolution in identifying microbial species and functional features but typically requires deep sequencing to detect unique species in complex communities. (Shapira et al., 2009) However, even shallow shotgun sequencing with as little as 0.5 million sequences per sample has been shown to be a powerful and costeffective alternative. (Hillmann et al., 2018) A significantly high correlation has been reported between 16S rRNA and shallow WMGS at both phylum and genus levels in our previous study.  Hence, we used shallow metagenomic data to identify disease-associated taxa and their functions. To achieve this, we performed a mixed linear regression model on metagenomics taxa (Supplementary Data S6) while adjusting for fixed effects of age, sex, BMI and random effects of house and recruitment site (Methods) in untreated MS versus HHCs, excluding the effect of treatment on gut microbiome. As shown in Figure 2A, compared to HHCs, 7 species were significantly reduced in untreated MS, including Fusicatenibacter saccharivorans, Blautia obeum and Faecalibacterium prausnitzii, whereas 16 species including Akkermansia muciniphila, Ruthenibacterium lactatiformans, Ruminococcus torques and Hugatella hathewayi were significantly increased in this group (FDR < 0.05). We observed a similar trend for these same species in untreated RRMS and progressive MS, although some did not reach significance likely due to the smaller sample size and relatively higher interindividual heterogeneity of these groups. Intriguingly, a larger decrease of F. saccharivorans and F. prausnitzii and a larger increase of R. lactatiformans, H. hathewayi and Eisenbergiella tayi were found in untreated PMS compared to untreated RRMS ( Figure 2B), suggesting the alteration of these species could be associated with disease progression.
To determine whether the microbial species identified in untreated RRMS or PMS were associated with disease severity, we tested the correlation between microbiota and the Multiple Sclerosis Severity score (MSSS), adjusting for age and BMI. Several species showed correlations with disease severity in untreated RRMS and PMS patients ( Figure 2C-D, Spearman's correlation, p < 0.05), consistent with a recent study. (Cox et al., 2021) Specifically, some Bacteroides species were correlated with lower MSSS in RRMS and short-chain fatty acid producers like Butyrivibrio, Clostridium and Ruminococcus were correlated with lower MSSS in PMS. Conversely, Collinsella aerofaciens, shown to increase disease severity in collagen-induced arthritis mice (Chen et al., 2016b) was associated with a higher MSSS in RRMS patients. Consistent with studies (Larsen, 2017) showing increased inflammatory properties of several Prevotella species (including P. buccalis, P. corporis, P. disiens, and P. copri) in chronic inflammation, we found these were associated with higher MSSS in progressive MS patients. However, the abundance of these species showed no significant difference in untreated RRMS (versus HHCs) or PMS (versus HHCs). Finally, Streptococcus thermophilus, Azospirillum sp. 47_25 and Rhodospirillum sp. UNK.MSG-17 were also correlated with MSSS albeit in different direction for RRMS (positive) and PMS (negative).

Functional alterations in the gut microbiome of untreated MS patients
To explore the functional potential of the MS gut metagenome, we applied the HUMAnN2 workflow to all individual samples. Generally, all microbes perform four core metabolic pathways, biosynthesis, degradation, energy metabolism and macromolecule modification ( Figure S3A). No significant differences in the functional potential of gut microbes were observed between MS and HHC. PCA analysis of the abundance of functional pathways also failed to identify significant changes between untreated MS patients (RRMS or PMS) and HHCs ( Figure S3B, PERMANOVA P > 0.05). However, when testing individual pathways, we found that phytate degradation I (PWY-4702), was significantly more represented in MS patients ( Figure 3A). Several species, including Akkermansia muciniphila, Escerichia coli and Cronobacter sakazakii, have the ability to degrade phytate via this pathway. We found the increase of both A. muciniphila ( Figure 2A) and two proteins, 4-phytase (Amuc_0145, B2ULU5, EC 3.1.3.26) and Inositol-1-monophosphatase (Amuc_1242, B2URI2, EC 3.1.3.25) from A. muciniphila involved in this pathway in untreated RRMS and progressive MS patients ( Figure 3B). As multiple (and sometime opposing) functional capabilities have been reported for A. muciniphila strains, we implemented the Metagenomic Intra-Species Diversity Analysis System (MIDAS)  to estimate strain-level genomic variation of A. muciniphila, including gene content and single nucleotide polymorphisms (SNPs), from our shallow shotgun metagenomic data. In total, 58 samples had sufficient sequencing coverage allowing us to identify single nucleotide polymorphisms (SNPs) and gene content from A. muciniphila ( Figure S4). To distinguish possible strains from these samples, we compared 2,913 A. muciniphila genes (presence/absence) and found four clusters of genes in two clusters of samples (cluster1 and cluster2). None of the sample clusters was significantly correlated with disease status, sex, treatment status or geographic site. The majority of A. muciniphila genes were shared across samples, (i.e. core genes), but some genes showed a distinct pattern. Functional analysis revealed that "Sulfite reductase [NADPH] hemoprotein beta-component (EC 1.8.1.2)", encoded by the cysI gene, was present in cluster1 but not cluster2. Intriguingly, Becken et al. reported the cysI gene as present in AmI phylogroup but missing in AmII and AmIV. (Becken et al., 2021) While additional studies are needed to establish their relevance to MS, we were able to identify the presence of at least two A. muciniphila strains with differences in sulfur metabolism.
As we observed some microbial species were significantly correlated with MS disease severity ( Figure 2C-D), we then tested which microbial pathways were involved. We found different pathways associated with disease severity in untreated RRMS and PMS ( Figure 3F). Particularly, abundances of "PWY-4981: L-proline biosynthesis II (from arginine)" was positively correlated with higher Multiple Sclerosis Severity score (MSSS) in untreated RRMS patients, and Collinsella aerofaciens is one the most dominant species involved in this pathway ( Figure S3C) which was found to be correlated with higher MSSS as well ( Figure 2C). Conversely, "PWY-5097: L-lysine biosynthesis VI" pathway, a pathway reduced in MS patients compared to HHCs ( Figure 3A), was related to a lower MSSS in untreated progressive MS dominantly contributing by species like Bacteroides species and F. prausnitzii ( Figure S3D).

Specific interacting microbial communities were enriched in MS
As the gut microbiome is an ecosystem in which microorganisms interact closely with each other, we next applied microbial co-abundance network inference to explore the interconnection pattern across microbes in untreated MS and their HHC. We computed species-species co-abundance relationships by using SparCC and controlled the significant correlations at FDR 0.05 level using 100 X permutation. In total, 1677 species were used for the analysis, resulting in 116,397 correlations across 1372 species in MS patients and 105,304 correlations across 1375 species in HHC (absolute Sparcc r ≥ 0.1, FDR < 0.05, Supplementary Data S7). To generate only highly correlated microbial networks, we filtered out correlations with r < 0.4 (based on the network centrality distribution) ( Figure S5A) and co-abundance subnetworks containing < 2 species. After correcting for age, sex and BMI differences between MS and HHC cohorts, we identified a network of 773 taxa with 5688 correlations in HHC dominated by 555 Firmicutes species and 196 Bacteroidetes species, and a network of 786 taxa with 6742 correlations in MS, dominantly by 549 Firmicutes species and 197 Bacteroidetes species ( Figure S5B). Notably, the majority of taxa (n = 702) and correlations (n = 4131) between MS and HHC microbial networks overlapped ( Figure  S5C), suggesting that the fundamental role of commensal microbes remains stable even under different biological conditions.
To test whether any part of the co-abundance microbial communities was specific for either group (MS or HHC), we performed a cohort-specific analysis (Methods) and found 215 correlations across 119 species specifically in untreated MS patients ( Figure 4A, mean r = 0.78, FDR < 0.05), and 195 correlations across 139 species specifically in HHCs ( Figure 4B, mean r = 0.783, FDR < 0.05). Species from the same phylum were clustered together in both MS and HHC networks, suggesting a cooperative role of these species in response to the environment. Remarkably, we observed different Firmicutes/Bacteroidetes (F/B) ratio for the MS-specific network (F/B=2.5) and HHC-specific network (F/B=1.03) (hypergeometric test p < 0.01). Interestingly, 45 unique species (largely dominated by Bacteroides and Prevotella species) composed the HHC network ( Figure  4C).
Next, we evaluated whether the differentially abundant species between MS and HHC contribute to these group -specific microbial networks. Surprisingly, among 21 significantly altered species (untreated MS versus HHC, Figure 2A), seven species, A. muciniphila CAG:154, Akkermansia sp. 54_46, Akkermansia sp. Phil8, Akkermansia sp. UNK.MGS-1, Tissierellia bacterium S7-1-4, Peptoniphilus grossensis and Porphyromonas bennonis, were identified in both MS and healthy specific networks, and only one species, Varibaculum cambriense, was found in the MS network ( Figure 4D). This suggests that co-abundance relationships and differential microbial abundance reflect orthogonal information. On the other hand, the group -specific species didn't show significant differences in abundance, but some have unique functions ( Figure 4E-F). For example, [Clostridum] innocuum and Salmonella enterica with unique functions mapped to drug resistance and pathogenicity were specific to the MS network, while Bacteroides vulgatus, Bacteroides thetaiotaomicron, Prevotella fusca and Prevotella denticola with unique functions mapped to glycan biosynthesis, were specific to the HHC network. Altogether, these results suggest that, compared to the alteration of single species, microbial co-abundance network analyses allowed us to identify a specific highly interacted community that may also be a contributor to health or disease status.

Impact of treatment on gut microbiota
We next evaluated how DMT may affect gut microbiome composition in RRMS patients receiving any of 6 commonly used treatments in our study, namely dimethyl fumarate (DMF, n=77), fingolimod (n=66), glatiramer acetate (GA, n=66), interferon (n=76), natalizumab (n=25) and anti-CD20 (rituximab and ocrelizumab, n=15). Overall, the microbial composition measured by βdiversity did not differ between treated and untreated RRMS patients (except for the interferon treated group). However, significant differences in β-diversity were observed when patients within each treatment group were compared to their corresponding HHC ( Figure 5A).
Due to the heterogeneity of treated and untreated RRMS patients recruited from multiple locations and unequal sample sizes of these groups, we mainly focused our analyses on gut microbiome by comparing untreated or treated RRMS groups to their HHCs ( Figure 5B-C). A direct comparison between untreated RRMS and treated RRMS was shown in Figure S6A-B. Intriguingly, the microbial changes in untreated RRMS patients (versus HHCs) were observed to show different changing patterns in treated RRMS (versus HHCs). Specifically, several taxa increased in untreated RRMS subjects showed no difference within DMTs groups, including Parabacteroides merdae CAG:48, A. muciniphila and other Akkermansia species. However, it is possible that the smaller n as a result of stratification may limit the statistical power to detect differences. Use of DMTs was also associated with changes in multiple taxa that were not significantly different between MS and HHC. For example, DMF, which is hydrolyzed into monomethyl fumarate (MMF) before exerting its therapeutic properties, specifically reduced Bacteroides stercoris, Clostridium and Eubacterium species, and fingolimod (a sphingosine-1-phosphate receptor modulator prevents the egress of activated lymphocytes from the lymphoid tissue) specifically reduced Bacteroides finegoldii CAG:203, Roseburia faecis and Blautia species. Interferon-ß treatment, thought to decrease proinflammatory cytokines and prevent the migration of activated T cells across the blood-brain barrier, was associated with dysbiosis of short-chain fatty acidproducing bacteria like Ruminococcus sp., Clostridium sp., Faecalibacterium prausnitzii, Roseburia inulinivorans and Roseburia intestinalis while also increased Parabacteroides distasonis, which have been shown to have multiple metabolic benefits in obesity.  Notably, Bacteroides uniformis was significantly increased by interferon treatment but reduced by glatiramer acetate and natalizumab therapy. This bacterium was reported to be associated with MS (Miller et al., 2015) but also with a protective role in obesity.(Lopez-Almela et al., 2021) GA exerted a modest impact on gut microbes compared to other DMTs. Lastly, infusion of natalizumab (an α4 integrin antagonist preventing leukocyte trafficking into the CNS) or anti-CD20 monoclonal antibody (causing B cell depletion) altered gut microbes significantly. Phascolarctobacterium sp. CAG:207 was increased while Prevotella species and Bifidobacterium longum were decreased in response to natalizumab. Reduction of Bacteroides finegoldii CAG:203 and Blautia sp. CAG:37 were observed in association with anti-CD20.
The alteration of microbial composition can lead to a change in their overall metabolic profile. Based on metagenomic sequencing, numerous metabolic pathways appeared to be altered under the DMTs, and many of them are included in the same high-class pathway ( Figure S6C). To better understand the outlined functions related to the treatment, we grouped the pathways by their hierarchical class and found that pathways related to lysine synthesis, sugar nucleotides and unsaturated fatty acids biosynthesis were decreased significantly in untreated RRMS but modulated differently by the various DMTs. Of interest, the increased cyclitols degradation pathways in untreated RRMS remained highly abundant even after treatment ( Figure 5C and Figure S6C). We also identified various metabolic pathways that were differentially modulated by specific therapies. For example, DMF use increased heme synthesis and enzyme cofactor biosynthesis pathways. In addition. DMF and interferon use was associated with an increase in Lornithine biosynthesis and carrier biosynthesis. Furthermore, GA use was associated with increased peptidoglycan biosynthesis and natalizumab with increased lipid biosynthesis, whereas a decrease of guanosine nucleotides degradation pathway was associated to Fingolimod treatment ( Figure 5C and Figure S6C). Altogether, disease modifying therapies showed significant and specific impact on gut microbiome both structurally and functionally, indicating the importance of stratifying microbiome analyses by treatment status.
To further investigate the mechanism of disease modifying therapies (DMT) in MS and their interactions with gut microbiota, we performed metabolomic profiling in 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. A panel of global metabolites and targeted short-chain fatty acids (SCFAs) in both feces and serum samples were measured using ultrahigh Performance Liquid Chromatography-Tandem Mass Spectroscopy (UPLC-MS/MS) (Supplementary Data S8).
We focused our analysis on the fecal metabolites involved in the reactions mediated by microbial proteins found by metagenomic sequencing and found 31 metabolites significantly different between untreated patients and controls, or in response to at least one MS drug ( Figure 6A). Consistent with their expected functions and origin, we found higher variability in fecal metabolites compared their corresponding serum levels (with the notable exception of increased serum fumarate in DMF-treated patients). We also identified significant changes in the levels of 8 SCFAs in either serum or stool for at least one group ( Figure 6B). Remarkably, the vast majority of changes in microbiota-derived fecal metabolites were towards lower levels among MS patients and even more significantly in response to DMTs (except for GA, Figure 6A). Higher levels of metabolic dysfunction have been reported to be associated with increased disability in MS. (Lazzarino et al., 2017;Villoslada et al., 2017) We found no difference of disease severity (measured by global MS severity score) among RRMS patients (treated or untreated) ( Figure 6C). This suggests the altered metabolites reported here are in response to the MS drugs, not the disease process. Interestingly, we found specific signatures of microbe-derived metabolites (stool) in RRMS patients in response to each treatment. The most notable changes in gut metabolites were induced in response to Fingolimod and IFN-β.
While Fingolimod is an oral drug, and changes to the gut microbiota might be expected, the profound metabolic signature of IFN-β (an injectable) was most intriguing. A functional analysis of the 23 IFN-β-associated metabolites, revealed a significant enrichment in pathways involving amino acid metabolism (e.g. "Arginine biosynthesis"), carbohydrate (i.e. "Citrate cycle"), nucleotide (i.e. "Purine"), and energy ("Nitrogen") metabolism ( Figure 6D). In contrast, GA exerted an almost null impact on stool-derived metabolites. These findings are in agreement with previous studies, in which only modest transcriptional changes were observed in PBMCs after treatment with GA compared to IFN-β.(De Jager et al., 2009;Ottoboni et al., 2012) Also of interest, these distinct metabolomic alterations were consistent with functional predictions derived from shotgun sequencing ( Figure 5B-C).
We next addressed the question of whether these changes in metabolite levels were linked to the relative levels of gut microbes as identified by metagenomics. We noted that pyruvate was significantly decreased in both feces and serum samples from RRMS patients treated with Fingolimod. Interestingly, this finding correlates well with the significant depletion of taxa containing the pathway "CARBOXYLATES-DEG" (which produces pyruvate) in Fingolimod treated patients ( Figure 5C, Figure S6). We also observed the concentration of fecal SCFA (such as acetate and propionate) was consistently lower in RRMS patients, regardless of treatment ( Figure 6B), consistent with our finding of the depletion of F. prausnitzii (a major SCFA-producing bacteria) in MS. A decreased amount of fecal SCFA has also been reported in RRMS and PPMS patients in other studies. (Takewaki et al., 2020;Zeng et al., 2019) Propionate supplementation in MS patients was associated with an increased Treg/Th17 ratio, leading to long-term clinical improvement. (Duscha et al., 2020) Interestingly, we found a significant increase in serum propionic acid (Acetic and Butyric acid also followed this same trend, without reaching statistical significance) in RRMS patients treated with IFN-β ( Figure 6E, Figure  6B). Since most SCFAs produced in the colonic lumen are actively transported to the lamina propria and further into the blood stream, (Olsson et al., 2021;Venegas et al., 2019) we hypothesized that IFN-β may increase the intestinal absorption of propionate, as part of its immunomodulatory effect. To address this hypothesis, we searched whether expression of the genes encoding for SCFA transporters MCT1 (SLC16A1) and SMCT1 (SLC5A8) (Miyauchi et al., 2004;Ritzhaupt et al., 1998) were upregulated by IFN-β.The Interferome database (Rusinova et al., 2013) reports an increase of SLC16A1 expression in human bronchial epithelial cells (no data is available for intestinal epithelial cells) treated with IFN-β ( Figure 6F), potentially supporting our hypothesis.

Diet and gut microbiome
Diet is thought to explain over 20% of microbial structural variations in humans, implying the potential for dietary strategies in disease management through gut microbiota modulation.  To explore the association of diet with host characteristics and the ability of diet to modulate gut bacteria in MS, we collected validated Block 2005 food frequency questionnaires (FFQs) (Block, 2005) from our participants. Among the 576 pairs of MS and HHC, 517 pairs (89.8%) finished the questionnaires, in which 37 nutrients were quantitated (Supplementary Data S2) through an external vendor (NutritionQuest). Recent epidemiologic studies of diet and health outcomes have also focused on the overall diet quality, (Guo et al., 2004) which can be measured by the Heathy Eating Index (HEI-2015), where a higher HEI-2015 score indicates greater diet quality (See Supplemental document).
We first evaluated whether diet differs significantly by disease course, household, recruitment site, and other variables of interest. As shown in Figure 7A, significant differences in diet were associated with BMI, participant household, recruitment site, education and age. Not surprisingly, a higher BMI correlated with a lower HEI-2015 score in both MS patients and healthy individuals ( Figure S7A), consistent with evidence that an imbalanced diet exerts a significant influence on weight. (Guo et al., 2004) We also observed that diet quality increased with age ( Figure S7A). There is considerable variation in dietary intakes across countries ( Figure 7B). In particular we found a significantly lower average HEI-2015 score in participants from Buenos Aires when compared to those in San Francisco, New York, Edinburgh and San Sabastian. While this may indeed indicate a lower health index, it is noteworthy that the FFQ is standardized for the US average participant, and diets in other parts of the world may not adjust properly to this standard. As gut microbial diversity differed among recruitment sites, we hypothesized that the diversity was influenced by diet. Indeed, we found that higher microbial diversity significantly correlated with a higher HEI-2015 score in both healthy and MS individuals ( Figure 7C, Pearson's correlation, p < 0.01). However, although participants from Buenos Aires had lower HEI-2015 scores, their microbial diversity remained high compared to other sites, whereas New York had higher HEI-2015 scores but comparatively lower microbial diversity ( Figure S2C). This may indicate that standardized questionnaires, even if validated, do not fully capture the wide range of dietary habits from iMSMS participants, but also suggests that the gut microbiota could be influenced by other factors, such as physiological activity, water and air, among other possible factors. Also, shifts in oral microbe composition need to be considered as studies have shown oralderived bacteria can colonize and persist in the intestines. (du Teil Espina et al., 2019;Hatton et al., 2018) Despite the large variance in dietary habits among participants, we identified a significantly higher diet similarity within household pairs compared to that of random pairs drawn from within the same city ( Figure S7B). The lowest diet similarity was found when random pairs of MS and HHC were assembled from different cities, consistent with our previous findings  and reflecting distinct dietary habits across cities and countries ( Figure 7B). Finally, we observed a significant correlation between education, nonsmoker (or former smoker) status, and female sex with a higher HEI-2015 score ( Figure S7C-E), also consistent with findings from previous studies. (Arabshahi et al., 2011;Thorpe et al., 2019) Although a more similar diet was shared among household participants, the HEI-2015 score of MS patients was significantly higher than those of healthy controls ( Figure 7D, paired T-test, p < 0.001). However, microbial taxa associated with MS status did not overlap with those associated with diet, thus likely not representing a confounder. Indeed, we specifically assessed which dietary components were consumed differently by MS and healthy participants and whether these differences were associated with species previously shown to be altered in MS. By comparing the ten components from the HEI-2015 (Supplementary Data S3), we found MS participants consumed more fruit, vegetables and unsaturated fatty acid when compared to HHCs ( Figure 7E), which contributed to their scores ( Figure 7D). A deeper analysis of the relationship between gut microbiota and diet led us to identify that Eubacterium eligens was highly correlated with a higher HEI-2015 score ( Figure 7F), and particularly correlated with intake of whole grains, fruit and vegetables (after adjusting for age, sex, BMI and recruitment site, Figure S7F), consistent with previous studies showing that E. eligens responded significantly to dietary fiber. (Chung et al., 2016) Faecalibacterium sp., Eubacterium sp. and Blautia sp. were also positively correlated with higher intake of whole grains. Increased Alistipes obesi abundance was also correlated with healthier diet ( Figure 7F). Interestingly, other studies found low Alistipes abundance in individuals with obesity (Thingholm et al., 2019) and was linked to higher meat ,(Garcia-Ribera et al., 2020) while it was identified as a predictor of successful weight loss in a two-year intervention (including healthier diet) in adult with obesity,  which suggests a potential beneficial role of Alistipes species in the context of metabolic health. Altogether, although diet does correlate with changes in the host microbiota, we were able to tease apart the effects of diet and disease in large part, thanks to our household paired design ( Figure S7F).
We further explored how taxa associated with untreated MS were related to diet. To exclude the impact of disease status or treatment, we stratified the analysis into healthy control, untreated MS and all sample groups. As expected, diet showed a modest effect on these MS-associated taxa after controlling for the environmental impact by household design in all three groups ( Figure 7G). Still, some disease-associated species were also related to diet. For example, Ruminococcus torques was enriched in MS, and showed a negative correlation with sodium intake, whereas no difference in sodium intake was found between MS versus HHCs. Faecalibacterium prausnitzii correlated positively with fruit (which MS patients consumed more), but the bacterium remained reduced in MS compared to healthy controls. These examples suggest that these species were more likely related to disease status than diet.
Phytate degradation I (PWY-4702) pathway was found to be overrepresented in MS patients ( Figure 2A). Phytate, a plant-based antioxidant compound, is a strong chelator of divalent minerals (e.g. calcium, magnesium, iron and zinc), which bacteria metabolize into myo-inositol, a compound with immunoregulatory properties, (Nerurkar et al., 2020) which was found at lower levels in MS sera and CSF. (Zahoor et al., 2021) Thus, we hypothesized that this bacterial pathway was activated: i) in response to increased dietary intake of divalent minerals by MS patients or; ii) as a compensatory mechanism to produce more myo-inositol. To test this hypothesis, we compared the dietary mineral intake between MS patients and their HHCs and found no significant difference (after adjusting for age, BMI and sex) ( Figure S8A). None of these minerals was significantly correlated with gut microbes or microbial pathways after multiple testing correction. Thus, we speculate the increased Phytate degradation pathway seen in MS may be a compensatory effect to restore normal myo-inositol levels.
Finally, we observed MS patients took more vitamin D supplementation than healthy controls ( Figure S8B), possibly in response to studies showing an association with reduced risk of developing MS and of disease activity in MS patients. (Munger et al., 2006;Runia et al., 2012) When assessing the impact of vitamin D usage on microbial composition, we were unable to find a correlation. We did find a trend towards a negative correlation with microbial α-diversity for both MS or HHCs samples, but without reaching significance ( Figure S8C). Similarly, β-diversity was not significantly influenced by vitamin D intake ( Figure S8D).

Discussion
In this paper we characterized the gut microbiome in 576 MS patients and genetically unrelated household healthy controls. Microbiome composition and function significantly differed across disease subtypes and responded differently to disease modifying treatments and were modestly associated with diet. We found that the microbial composition was to a lesser extent, associated with factors such as geographic location, age, sex and BMI. The influence of other confounding factors was reduced by our paired household design, thereby potentially enhancing power to identify MS-associated microbial features. In addition to confirming and extending previous reports, (Berer et al., 2017;Cekanaviciute et al., 2017;Chen et al., 2016a;Cox et al., 2021;Jangi et al., 2016) this work provides a large reference dataset that can be used to understand microbial variation across individuals with MS, disease subtypes and in response to different therapeutic interventions.
When studying abnormalities of microbial composition, a common approach is to evaluate changes in -diversity and -diversity. Consistent with earlier studies, we found no difference of diversity between MS patients and healthy individuals (Berer et al., 2017;Cekanaviciute et al., 2017;Jangi et al., 2016) (either when stratified as untreated MS versus HHCs or as treated MS versus HHCs). However, in contrast to previous studies, we observed a significant difference of -diversity in disease status (regardless of treatment status) compared to healthy controls. Interestingly, there was no difference in -diversity between untreated MS and treated MS, potentially indicating that disease status exerts a stronger effect on gut microbiome than does treatment. (Cox et al., 2021) Overall, our findings revealed a robust alteration of gut microbial composition related to the disease and therapy.
In line with our previous report,  recruitment site explained the largest proportion of variance in microbial composition, followed by other confounders such as age, sex and BMI. These results support the concept that large, multicenter studies provide greater statistical power while retaining the ability to control for the main confounders. (Abeles et al., 2016;Goodrich et al., 2014;Lax et al., 2014;Song et al., 2013) While an increase in A. muciniphila has also been reported in previous studies, (Berer et al., 2017;Cekanaviciute et al., 2017;Cox et al., 2021;Probstel et al., 2020) interpretation of its specific role remains controversial. A. muciniphila is a mucin-degrading bacteria shown to exert proinflammatory effects on T cells in vitro  and to exacerbate inflammation during infection. (Ganesh et al., 2013) Interestingly, peptides from A. muciniphila have been recently shown to stimulate human myelin autoreactive CD4+ T cell clones, thus suggesting molecular mimicry is a potential mechanism for MS pathogenesis . However, A. muciniphila has also been proposed as a contributor to maintaining gut health, improving glucose homeostasis, increasing gut mucin integrity and enhancing effect of checkpoint inhibitor immunotherapy. (Cani and de Vos, 2017;Liu et al., 2019;Routy et al., 2018) Different functional capabilities across A. muciniphila strains that may affect how these bacteria interact with the host. (Becken et al., 2021;Karcher et al., 2021;Kirmiz et al., 2020) At least two A. muciniphila strains were identified in our samples with differences in their functions such as sulfur metabolism, but none of them was enriched in MS in our dataset. Through pathway analysis we found "phytate degradation I" (PWY-4702) (a cyclitols degradation pathway), mainly driven by A. muciniphila, was significantly increased in untreated MS patients. This pathway converts phytate into Myoinositol. Phytate is a strong chelator of divalent minerals such as calcium, magnesium, iron and zinc. Previous studies suggested that high levels of iron and zinc could play a role in MS activity and progression, (Ferreira et al., 2017;Hametner et al., 2013;Sanna et al., 2018) whereas calcium and magnesium could exert neuroprotective capacities. (Enders et al., 2020;Goldberg et al., 1986) Dietary mineral intake was no different between MS and healthy controls, but it is still possible that bacterial pathways (such as Phytate degradation) modulate the bioavailability of these minerals, thus contributing to disease pathogenesis. Myo-inositol, a simple carbohydrate produced in the body and available in foods such as fruits and cereals, is involved in lipid signaling, osmolarity, glucose, and insulin metabolism (Gonzalez-Uarquin et al., 2020) and utilized as dietary supplementation in different pathological conditions, including diabetes and metabolic disorders. (Pintaudi et al., 2016;Shokrpour et al., 2019) Interestingly, a very early study showed that patients with MS appeared to metabolize myo-inositol abnormally, (Holm, 1978) and administered myo-inositol was shown to have a positive effect on evoked potential responses in MS (n =9) and controls (n=9). (Young et al., 1986) The role of Akkermansia in myo-inositol metabolism needs to be further elucidated.
Ruminococcus torques is another potent mucus degrader and may decrease gut barrier integrity. (Cani, 2014;Rajilic-Stojanovic and de Vos, 2014) A recent study showed that R. torques was associated with an enhanced MRI T2 signal in multiple motor brain areas and exacerbated disease in an animal model of amyotrophic lateral sclerosis (ALS). (Blacher et al., 2019) Ruthenibacterium lactatiformans, a lactate-producing species, was previously associated with an increased EDSS and decreased lower extremity motor function in RRMS and progressive MS. (Cox et al., 2021) Overall, 7 species were significantly reduced in untreated MS. Faecalibacterium prausnitzii, one of the main butyrate producers found in the intestine, has anti-inflammatory properties that were partly associated with secreted metabolites that block NF-κB activation, IL-8 production and upregulate regulatory T cell production.  It can also attenuate the severity of inflammation through release of metabolites that enhance intestinal barrier function. (Carlsson et al., 2013;Martin et al., 2015) The pyruvate-producing carboxylates metabolism pathways, contributed by F. prausnitzii, were found to be significantly reduced in untreated MS patients. Altogether, we found a depletion of potentially beneficial bacteria in untreated MS patients compared to healthy controls, which in turn disturbed key metabolic pathways that might be expected to worsen the inflammation of MS. These findings could lead to the development of "designer probiotics" that can restore the healthy composition and function of the gut microbiome.
We next tested whether these altered bacteria also associated with disease severity, and found that only Streptococcus thermophilus, Azospirillum sp. 47_25 and Rhodospirillum sp. UNK.MSG-17 were. However, correlations were positive for RRMS and negative for PMS patients. This implies that the change of gut microbial community may be linked to the onset of disease and stabilized during the disease course, a hypothesis which requires further investigation by longitudinal studies. Several other species were found to be correlated exclusively with MS severity (e.g. not with disease status). For example, Butyrivibrio, Clostridium and Ruminococcus species, which are short-chain fatty acid (SCFA) producers, correlated with lower MS severity in PMS. It's well known that SCFAs play a critical role in immunoregulation with well-characterized anti-inflammatory effects on both epithelium and peripheral immune cells. This implies potentially beneficial effects of these bacteria by producing anti-inflammatory metabolites. On the other hand, Collinsella aerofaciens, a species showed to increase disease severity in collagen-induced arthritis mice, (Chen et al., 2016b) was associated with a higher MSSS in RRMS patients probably via the pathway "PWY-4981: L-proline biosynthesis II (from arginine)". Prevotella species such as P. buccalis, P. corporis, P. disiens, P. copri were associated with higher MSSS. Although Prevotella species have been associated with health-beneficial properties, several studies have shown associations with autoimmune diseases, insulin resistance and diabetes, and gut inflammation. (Leite et al., 2017;Pedersen et al., 2016;Scher et al., 2013) Intriguingly, we found pathway "PWY-5097: L-lysine biosynthesis VI", a decreased pathway in MS versus HHCs, was associated with a lower disease severity. Several commensal bacteria participants in this pathway, including Bacteroides, Faecalibacterium and Eubacterium species. L-lysine has been shown to have anti-inflammatory in rat with chronic lung injury, (Zhang et al., 2019) and may play a neuroprotective role in intracerebral hemorrhage injury, (Cheng et al., 2020) thus suggesting a potential usage of L-lysine to suppress the disease progression. Based on these observations, we speculate that the role of gut bacteria in disease progression/severity is multi-faceted and individual-dependent.
Although structural and functional changes of gut microbiota in untreated MS patients have been identified in multiple studies, very little is known about how DMTs modulate these microbial communities. In this large dataset, we were able to study the effects of 6 DMTs including oral (Fingolimod, DMF), subcutaneous (GA, IFN) and infusion (Natalizumab, B-cell depleting) therapies. The MS associated alterations in gut microbiota may be modulated by the therapies as reported in other studies. (Jangi et al., 2016;Sand et al., 2019;Storm-Larsen et al., 2019) This is supported by our finding that several species showing differential abundance between RRMS and HHCs (Firmicutes bacterium CAG:65, Sutterella wadsworthensis, Parabacteroides merdae CAG:48, Akkermansia muciniphila and others), showed no difference in treated groups versus HHCs. In addition, DMTs resulted in a decrease in the relative abundance of specific taxa that are not MS-associated, potentially by their innate antimicrobial properties. (Maier et al., 2018;Storm-Larsen et al., 2019) Compared to untreated RRMS, several common gut microbes including Bacteroides, Blautia and Clostridium species were significantly reduced in response to oral medications and species like Faecalibacterium prausnitzii, Dialister invisus CAG:218 and Roseburia intestinalis were reduced in individual receiving injectables. Furthermore, infused therapies resulted in a decrease of Bifidobacterium adolescentis, which was shown to promote Th17 cell accumulation and exacerbated autoimmune arthritis in a mouse model, arguing for its pathological relevance. (Tan et al., 2016) On the other hand, we found several species that were increased by DMTs, in particular Ruthenibacterium lactatiformans and Ruminococcus torques (with Fingolimod), Eubacterium hallii (with GA) and Bacteroides uniformis (with interferon). Intriguingly, sequence-based analysis suggested the oral drug Fingolimod would induce the most metabolic changes compared to other medications, a finding validated by metabolomic analysis. Specifically, the depletion of microbial "CARBOXYLATES-DEG" pathways (which produces pyruvate) may explain the low level of pyruvate observed in feces and serum samples from RRMS patients treated with Fingolimod, and the depletion of F. prausnitzii (a major SCFA producing bacteria) could account for the lower levels of acetate and propionate found in MS. We also found several microbe-derived fecal metabolites were remarkably lower in treated RRMS patients, implying a particularly important effect of these medications, likely through direct interactions with gut microbiota. Of interest, a significant increase of serum propionic acid was found in RRMS patients treated with interferon. Propionate supplementation in MS patients has been associated with an increased Treg/Th17 ratio, leading to long-term clinical improvement. (Duscha et al., 2020) Based on our findings, we propose the increased absorption of microbially-derived propionate via upregulation of the SCFA transporter MCT1 (SLC16A1) as contributing mechanism of action for IFN-β. Our results provide compelling evidence that DMTs have considerable effects on gut microbiota, not only compositionally but functionally, that may highlight therapeutic mechanisms requiring further investigation. Additional larger and longitudinal follow-up studies will help to evaluate these effects more precisely.
We demonstrated that diet quality (as measured by HEI-2015), was significantly associated with BMI, age and geographic location. Furthermore, a higher education, nonsmoker or former smoker status, and being female were also associated with the higher HEI-2015, consistent with previous studies. (Arabshahi et al., 2011;Thorpe et al., 2019) A healthier diet associates with higher microbial diversity, but diet may not the only factor at play. Some bacteria remained unaffected by dietary change depending on host phenotype and the preexisting microbiota composition. (Flandroy et al., 2018) In addition, local environment (i.e. air, soil and water) could also influence diversity of the gut microbiota by horizontal transmission of environmental microbes. (Tasnim et al., 2017) Due to shared environment and a similar diet taken by the household MS and healthy individuals, we ensured that the household design can not only effectively control the impact of diet, but also the influence of environmental microbes on gut microbiota. Even when sharing a diet with healthy individuals, MS patients tend to eat heathier by taking more fruit, vegetables and unsaturated fatty acids. Vitamin D deficiency has long been associated with MS, and higher levels of vitamin D were associated with reduced clinical activity in established MS. (Munger and Ascherio, 2011) Unsurprisingly, we observed that MS patients took more vitamin D, but showed no significant influence on gut microbiome composition.

Limitations of the study
Shotgun metagenomics sequencing was limited to ~500,000 reads per sample. While this coverage is adequate to classify bacterial communities with higher resolution that 16S RNA gene sequencing, and to provide some insight into the metabolic potential of the communities, higher sequencing depth will be needed to resolve most strains, clades and DNA polymorphisms. We cannot exclude power limitations due to stratification by treatment. As a consequence of the paired household design, the majority of the pairs are spouses, thus leading to an uneven sex distribution of MS (69.4% of the MS participants are female, in keeping with the expected demographics for MS (Langer-Gould et al., 2013)). However, our model adjusted for the effect of sex on gut microbiome.
In summary, this is a large, multi-center gut microbiome study in MS patients and HHC. The findings strongly support the presence of specific gut microbiome associations both with MS disease course and progression, and functional changes in response to treatment. The origin and biological relevance of these associations remain to be elucidated. Nevertheless, our study supports the possibility that microbial manipulation and dietary intervention could be used as preventive and therapeutic strategies in MS.

Consortia
The members of iMSMS are: Xiaoyuan Zhou, Ryan Baumann, Xiaohui Gao, Myra Mendoza, Sneha Singh, Ilana Katz Sand, Zongqi Xia, Laura M. Cox, Tanuja Chitnis, Hongsup Yoon, Laura Moles, Stacy J. Caillier, Adam Santaniello, Gail Ackermann, Adil Harroud, Robin Lincoln, Refujia Gomez, Antonio González Peña, Elise Digga, Daniel Joseph Hakim, Yoshiki Vazquez-Baeza, Karthik Soman, Greg Humphrey, Mauricio Farez, Lisa Ann Gerdes, Jorge R. Oksenberg, Scott S. Zamvil, Siddharthan Chandran, Peter Connick, David Otaegui, Tamara Castillo-Triviño, Stephen L. Hauser, Jeffrey M. Gelfand, Howard L. Weiner, Reinhard Hohlfeld, Hartmut Wekerle, Jennifer Graves, Amit Bar-Or, Bruce A.C. Cree, Jorge Correale, Rob Knight, Sergio E. Baranzini (See Document S1 for consortium member affiliations).    . Each node indicates one species and color indicates the phylum classification. Each edge represents one species-species co-abundance relationship. (C) Overlapped counts of species and co-abundances in untreated MS specific and HHC specific networks. (D) Differential species in untreated MS versus HHC were overlapped with cohort specific species. (E-F) Functional pathways unique to the species highlighted in untreated MS specific network (E) and HHC specific one (F). Line size indicates the betweenness centrality of a species in the cohort specific co-abundance network.      Blue squares indicate that a gene is present, and yellow squares indicate that a gene was absent. Function was annotated for the genes missing in cluster2.

Figure S5. Microbial co-abundance in untreated MS and household healthy controls. (A)
Normalized centrality distribution of species-species co-abundance networks constructed in untreated MS and household healthy controls (HHCs) by setting correlation cutoff from 0.1 to 1 (SparCC, FDR < 0.05). (B) Microbial co-abundance network built in healthy (left, 773 species, 5688 co-abundances) and untreated MS individuals (right, 786 species, 6742 coabundances) by SparCC |r| ≥ 0.4 and FDR < 0.05 after adjusting age, sex and BMI. Each node indicates a species and color indicates the family classification. Each edge represents one species-species co-abundance relationship and labeled in green for positive correlation, red for negative correlation. (C) Overlapped counts of species and coabundances in untreated MS and HHC microbial networks (SparCC |r| ≥ 0.4, FDR < 0.05 adjusted for age, sex and BMI).

Supplement Document 1
Supplemental Document S1

Recruitment.
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 3 collection vials -a Q-tip (Q, dry), a snap frozen vial (S, wet), and a vial filled with Luria-Bertani broth (LB) and 30% glycerol. Participants were instructed to freeze the samples for at least 12 hours and ship them frozen with the ice pack included in the kit. Samples were returned to each site via overnight shipping in a thermal envelope. All participants were required to complete a clinical survey to report the disease status and treatment, and a subject survey to report demographic, medication, lifestyle and physiology factors. Clinical outcomes included the Expanded Disability Status Scale (EDSS), (Kurtzke, 1983) and the Multiple Sclerosis Functional Composite (MSFC). (Fischer et al., 1999) All data were collected and stored through secure REDCap TM questionnaires.

Sample preparation for sequencing.
To test whether the DNA processing method changes microbial composition, we extracted DNAs from the same 20 samples using both QIAcube and epMotion platforms. A subset of 40 samples prepared in Cohort 1 were re-sequenced in Cohort 2 to test the impact of sequencing runs on microbial composition.

Microbial diversity.
Since the MS and control subjects within household are often of different sex, the random comparisons between households utilized only cross-sexual comparisons to control for the effect of gender.

Shallow whole metagenome shotgun sequencing (WMGS) and data processing.
To deal with sparse microbial data in the analysis, we focused on species present in at least 5% of samples, covering at least 10 total reads. This provided a list of 1,677 species for use in the statistical analysis. The relative abundances of gene families were characterized from UniProt Reference Clusters (UniRef90) using HUMAnN2 (V2.8.2), (Franzosa et al., 2018) which were further mapped to microbial pathways and high-classes based on pathway hierarchy from the MetaCyc metabolic pathway database. (Caspi et al., 2016;Caspi et al., 2018) 490 pathways present in at least 5% of samples were retained for statistical analysis. Microbial gene families present in more than 5% samples were used to link with select fecal metabolites. The phylogenetic diversity of Akkermansia muciniphila was measured using Metagenomic Intraspecies Diversity Analysis System (MIDAS)  with its default parameters.

Diet analysis
Healthy Eating Index-2015(HEI-2015) was used for evaluation of the diet quality and calculated by NutritionQuest. The HEI-2015 adequate dietary components include 'total fruit', 'whole fruit', 'total vegetables', 'greens and beans', 'whole grains', 'dairy', 'total protein', 'seafood and plant proteins', and 'fatty acids', which are recommended to be high

SUPPLEMENTAL DISCUSSION
To identify disease -associated microbial features, a large microbiome study with participants recruited from multiple sites will increases statistical power to detect disease associated features, but also can introduce a greater number of confounding factors such as geography, host lifestyle, diet, and medication use that could shape the gut microbiome. (Gupta et al., 2017; Using shotgun metagenomics and adjusting for household, age, sex and BMI, we identified 16 species that were significantly increased in untreated MS versus HHCs. The same trend was observed for untreated RRMS (versus HHCs) and untreated progressive MS (versus HHCs), but some of these did not reach significance. This could be explained by the reduced power due to stratification. Specifically, Akkermansia muciniphila, Ruminococcus torques, Ruthenibacterium lactatiformans and Hugatella hathewayi were increased in MS. Hugatella hathewayi, (previously known as Clostridium hathewayi), was associated with several disorders, including acute appendicitis. (Randazzo et al., 2015) Interestingly, we found that H. hathewayi and another species Staphylococcus sp. CAG:324 increased even further in PMS compared to RRMS. By contrast, Dialister invisus was decreased in PMS but increased in untreated RRMS; This bacteria has been isolated from both human oral cavity and gut, and is considered a significant human pathogen. (Rocas and Siqueira, 2006) Several other studies have reported a decrease of F. prausnitzii in MS patients (Miyake et al., 2015;Swidsinski et al., 2017) and also in different intestinal disorders,  suggesting this organism may have beneficial properties. F. saccharivorans belongs to Clostridium subcluster XIVa and plays a critical role in immune system homeostasis by inducing regulatory T cells (Tregs) through the production of butyrate and other short chain fatty acids. (Murakami et al., 2018) Interestingly, a depletion of Clostridium subcluster XIVa has been reported in MS patients. (Takeshita et al., 2016) Also, low abundance of F. saccharivorans was found in active ulcerative colitis patients, while the administration of F. saccharivorans ameliorated oxazolone-induced colitis in a mouse model via suppression of IL-13 production. (Takeshita et al., 2016) Blautia obeum, an acetate-producing species, can inhibit the proliferation of Clostridium perfringens, Vibrio cholerae and vancomycin-resistant enterococci, also suggesting a possible therapeutic role.  Microbial co-abundance networks were computed to better understand how MS-associated microbes could drive changes in microbial communities. We observed a substantial alteration of interconnected species with an increased Firmicutes to Bacteroidetes (F/B) ratio in the untreated MS specific network. An increased F/R ratio is considered a hallmark of obesity and has also been associated with several pathological conditions (Magne et al., 2020) and we hypothesize the overrepresentation of Firmicutes in the MS-associated microbial network may contribute to pathogenesis. Surprisingly, few species individually associated with MS were part of the coabundant microbial networks, implying that microbial-induced pathogenesis may not be driven by individual differences level but instead by shifts community structure.
Diet is a key modifiable factor influencing the composition of the gut microbiota, indicating the potential for therapeutic dietary strategies to manipulate the microbial dysbiosis related to diseases. A number of human and animal feeding studies have provided insight into the effect of specific foods or dietary components on the gut microbial composition. (David et al., 2014;Li et al., 2009) We comprehensively explored the ability of the host diet to modulate not only gut bacteria in healthy individuals but also differential bacteria in MS patients.
An increased risk of MS was found to be associated with high energy and saturated fat intake, while food components, including fruit, vegetables, whole grains and unsaturated fats, appeared to exert a protective effect.(Katz Sand, 2018) Many MS patients subscribe to a special diet or dietary supplements. Intriguingly, we found that a healthier diet correlated with an increase in beneficial bacteria like Eubacterium and Faecalibacterium species. We also found that the differential species between MS and HHCs were unlikely to be shaped by diet. This may be explained by the household design that increases the power to identify true disease associated microbes by controlling the impact of environment.