DOI: 10.7272/q60c4t26
title: Gut microbiome of multiple sclerosis patients and paired household healthy controls reveal associations with disease risk and course

description: 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.
description: 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<sup>TM</sup>).  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 <em>SpiecEasi</em> 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 <em>cor2por</em> function from R package
 “<em>corpcor</em>”. 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). <strong>Global metabolomic profiling
 </strong> 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. <strong>Targeted short-chain fatty acid
 profiling</strong> 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<strong> </strong>nutrient<strong>
 </strong>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).