A Metagenome-Wide Association Study of Gut Microbiome in Patients With Multiple Sclerosis Revealed Novel Disease Pathology

Abstract

While microbiome plays key roles in the etiology of multiple sclerosis (MS), its mechanism remains elusive. Here, we conducted a comprehensive metagenome-wide association study (MWAS) of the relapsing-remitting MS gut microbiome (n_case = 26, n_control = 77) in the Japanese population, by using whole-genome shotgun sequencing. Our MWAS consisted of three major bioinformatic analytic pipelines (phylogenetic analysis, functional gene analysis, and pathway analysis). Phylogenetic case-control association tests showed discrepancies of eight clades, most of which were related to the immune system (false discovery rate [FDR] < 0.10; e.g., Erysipelatoclostridium_sp. and Gemella morbillorum). Gene association tests found an increased abundance of one putative dehydrogenase gene (Clo1100_2356) and one ABC transporter related gene (Mahau_1952) in the MS metagenome compared with controls (FDR < 0.1). Molecular pathway analysis of the microbiome gene case-control comparisons identified enrichment of multiple Gene Ontology terms, with the most significant enrichment on cell outer membrane (P = 1.5 × 10^-7). Interaction between the metagenome and host genome was identified by comparing biological pathway enrichment between the MS MWAS and the MS genome-wide association study (GWAS) results (i.e., MWAS-GWAS interaction). No apparent discrepancies in alpha or beta diversities of metagenome were found between MS cases and controls. Our shotgun sequencing-based MWAS highlights novel characteristics of the MS gut microbiome and its interaction with host genome, which contributes to our understanding of the microbiome's role in MS pathophysiology.

Keywords: dysbiosis; genome-wide association study; gut microbiome; metagenome shotgun sequencing; multiple sclerosis.

Figure 1

MWAS results of MS case-control phylogenetic association tests. (A) A quantile-quantile plot of the MWAS P-values of the clades. The x-axis indicates empirically estimated median -log₁₀ P-values. The y-axis indicates observed -log₁₀ P-values. The diagonal gray line represents y = x, which corresponds to the null hypothesis. The horizontal red line indicates the empirical Bonferroni-corrected threshold (α = 0.05), and the brown line indicates the empirically estimated (FDR-q = 0.1). Clades with FDR-q < 0.1 are plotted as red dots, and other clades as black dots. (B) A volcano plot. The x-axis indicates effect sizes of generalized linear model. The y-axis, horizontal lines, and dot colors are the same as in panel (A). (C) Phylogenetic tree. Levels L2–L7 are from the inside layer to the outside layer. The size and color of dots represent relative abundance and effect sizes, respectively. The five clades with suggestive case-control associations (FDR-q < 0.1) are outlined in red. FDR, false discovery rate; MWAS, metagenome-wide association study; MS, multiple sclerosis.

Figure 2

MWAS results of MS case-control gene association tests. (A) A quantile-quantile plot (left) and a volcano plot (right) of the MWAS P-values of genes based on the UniRef90 protein database. (B) A quantile-quantile plot (left) and a volcano plot (right) of genes based on the KEGG gene database. In the quantile-quantile plots, the x-axis indicates empirically estimated median -log₁₀ P-values. In the volcano plot, the x-axis indicates beta of generalized linear model as effect sizes. The y-axis in both plots indicates observed -log₁₀ P-values. The diagonal gray line represents y = x, which corresponds to the null hypothesis. The horizontal red line indicates the empirical Bonferroni-corrected threshold (α = 0.05), and the brown line indicates the empirically estimated (FDR-q = 0.1). Genes with FDR-q < 0.1 are plotted as red dots, and other genes as black dots.

Figure 3

MWAS results of MS case-control pathway association tests. (A) A quantile-quantile plot of the MWAS P-values of enrichment analyses based on GO terms. GO terms with P-values less than Bonferroni thresholds are plotted as red dots, and the other clades as black dots. (B) A quantile-quantile plot of the MWAS P-values of enrichment analyses based on KEGG pathways. (C) Comparison of P-values of GO enrichment analyses between the MS MWAS and GWAS data. The x-axis indicates the P-values of the GWAS. The y-axis indicates the P-values of the MWAS. The horizontal and vertical black lines indicate P-value of 0.05. The overlap of the GO enrichment was evaluated by classifying the GO terms based on the significance threshold of P < 0.05 or P ≥ 0.05 and using Fisher’s exact test. GO, Gene Ontology; GWAS, genome-wide association study; KEGG, Kyoto Encyclopedia of Genes and Genomes; MWAS, metagenome-wide association study; MS, Multiple sclerosis.

Figure 4

MS case-control comparison of microbial diversities. (A) Alpha diversities of the phylogenetic relative abundance data for six levels. Welch’s t-test of Shannon index between MS cases and controls showed no significant difference at any level. (B) Alpha diversities of the gene abundance data of the UniRef90 protein and KEGG gene databases. No significant case-control difference was found. (C) Beta diversities of phylogenetic relative abundance data at six levels. PERMANOVA based on Bray-Curtis dissimilarities found no significant differences among levels for either sequencing group with Bonferroni correction. (D) Beta diversities of the gene abundance of the UniRef90 protein database. No significant case-control difference was found. KEGG, Kyoto Encyclopedia of Genes and Genomes; NMDS, non-metric multidimensional scaling; PERMANOVA, permutational multivariate analysis of variance; MS, Multiple sclerosis.