Breakdown of multiple sclerosis genetics to identify an integrated disease network and potential variant mechanisms

Abstract

Genetics of multiple sclerosis (MS) are highly polygenic with few insights into mechanistic associations with pathology. In this study, we assessed MS genetics through linkage disequilibrium and missense variant interpretation to yield a MS gene network. This network of 96 genes was taken through pathway analysis, tissue expression profiles, single cell expression segregation, expression quantitative trait loci (eQTLs), genome annotations, transcription factor (TF) binding profiles, structural genome looping, and overlap with additional associated genetic traits. This work revealed immune system dysfunction, nerve cell myelination, energetic control, transcriptional regulation, and variants that overlap multiple autoimmune disorders. Tissue-specific expression and eQTLs of MS genes implicate multiple immune cell types including macrophages, neutrophils, and T cells, while the genes in neural cell types enrich for oligodendrocyte and myelin sheath biology. There are eQTLs in linkage with lead MS variants in 25 genes including the multitissue eQTL, rs9271640, for HLA-DRB1/DRB5. Using multiple functional genomic databases, we identified noncoding variants that disrupt TF binding for GABPA, CTCF, EGR1, YY1, SPI1, CLOCK, ARNTL, BACH1, and GFI1. Overall, this paper suggests multiple genetic mechanisms for MS associated variants while highlighting the importance of a systems biology and network approach when elucidating intersections of the immune and nervous system.

Keywords: GWAS; data integration; eQTL; multiple sclerosis; omics.

Fig. 1.

Workflow used to identify and study genes/variants involved in multiple sclerosis.

Fig. 2.

STRING Network and Gene Ontology (GO) Enrichment of multiple sclerosis genes. A network was created using the genes in Table 1. GO terms significantly enriched are shown with number of genes mapped over total with term across the genome and the false discovery rate (FDR) listed for each.

Fig. 3.

Missense variants from genome-wide association study (GWAS) and ClinVar for multiple sclerosis in CD6, EVI5, DKKL1, NR1H3, and HNRNPA1. A: deep codon evolutionary analysis performed on each gene (CD6, EVI5, DKKL1, NR1H3, HNRNPA1) with the number of species shown in parentheses and the annotated domains shown below. B: ranking of all variants from ClinVar (benign = yellow, pathogenic = red, variants of uncertain significance VUS = magenta) and gnomAD (gray). C: combined impact using gnomAD allele counts by variant scores for only variants at sites under conservation or selection. D–G: CD6, relative predicted impact of all gnomAD variants for CD6 on the top right identifying T217M as the highest impact (D). E: conservation of amino acids near site 217, which is in the middle of a potential phosphorylation site. F: model of CD6 with ConSurf coloring of conservation (yellow = conserved hydrophobic, blue = conserved basic, red = conserved acidic, green = conserved polar, gray = not conserved) embedded within a lipid membrane (cyan) and a zoom in view of site 217 on the right. G: the other linkage disequilibrium (LD) missense single nucleotide polymorphism (SNP) for CD6 showing no to weak conservation near amino acid 606. H: EVI5, conservation of amino acids around site 623 that are found throughout evolution as multiple amino acids. I: DKKL1, conservation around amino acid 109 that has no conservation. J: NR1H3, ClinVar variant at site 415, which is found in a highly conserved region of the protein. K: HNRNPA1, assessment of all ClinVar variants for the gene using multiple tools and our conservation analysis.

Fig. 4.

Human expression of genes within the multiple sclerosis list. A: a z-score metric of tissues within human FANTOM database addressing the number of multiple sclerosis (MS)-associated genes within that tissue with a z-score >2 (y-axis, i.e., with expression levels at least two standard deviations above the mean) and the average z-score for the MS genes in the tissues (x-axis). Genes mapped to the top right are labeled as they are likely tissues with highest expression of MS genes. B: expression z-scores of candidate genes in neutrophils (y-axis) and medial temporal gyrus (x-axis) segregating genes of the immune system and neurons. C: heat map for average expression of each MS-associated gene (Table 1) throughout 20 single cell organ data sets. Boxed in black are genes primarily associated with non-myeloid brain cells and those in orange with the immune system. D: gene clustering correlations in the thymus (y-axis) or spleen (x-axis) for cells that do or do not express MS genes. E: genes enriched in expression (y-axis) or number of cells (x-axis) similar to MS genes in the brain nonmyeloid data set.

Fig. 5.

Expression quantitative trait locus (eQTL) analysis of multiple sclerosis genome-wide association study (GWAS) linkage disequilibrium (LD) single nucleotide polymorphisms (SNPs) to nominate genes and epigenetic mechanisms. A: the % of MS genes in the list with eQTLs in multiple human tissues (gray), and the normalized value relative to the total number of eQTLs (including non-MS-associated genes) mapped within each tissue (red). B: genes with the most tissues having eQTLs. C: effect size (x-axis) and P value (y-axis) of eQTLs in subcutaneous adipose, highlighting the genes with the largest effect size. D: view of functional SNPs near the HLA-DRB1/5 regions. Shown on the top is Cohesin Hi-ChIP looping data of this site interacting with distal regions. E: ChromHMM models of several tissues showing the region to be a potential enhancer. Shown below are LD SNP correlations.

Fig. 6.

RegulomeDB of functional noncoding SNPs in the rs10492972 region. Starting from the top is the looping data for the rs61784580 promoter SNP or the intronic lead SNP rs10492972 followed by multiple data sets of the human genome browser, 18 state ChromHMM genome annotation for the region, CEU LD SNP correlation, details of transcription factor binding sites near both sites, and the eQTLs for KIF1B from the LD block all showing similar effect size.

Fig. 7.

Multiple sclerosis (MS) linkage disequilibrium (LD) single nucleotide polymorphisms (SNPs) overlap with other genome-wide association study (GWAS) traits. A: number of SNPs mapped per trait that overlap with any of the MS-associated LD block variants. Traits are ranked according to occurrence. B: SNPs that have the most traits associated with them that are also within the MS LD blocks. Top four are detailed in C, with two of the SNPs in LD with 0.93R² correlation. C: breakdown of three top SNPs for traits. For each, the far left shows the traits with the 95% confidence interval. Next is the analysis of functional LD SNPs through RegulomeDB showing the R2 value for each SNP to the lead SNP and the RegulomeDB score for the variant. After this is the list of eQTLs for the top RegulomeDB SNP. The far right shows expression data for the gene with the most significant eQTL from either FANTOM (immune cells) or Mouse Cell Atlas (brain).

Fig. 8.

Working model of pathways associated with multiple sclerosis (MS) etiology based on our genetic systems data integration.