Locus for severity implicates CNS resilience in progression of multiple sclerosis

Abstract

Multiple sclerosis (MS) is an autoimmune disease of the central nervous system (CNS) that results in significant neurodegeneration in the majority of those affected and is a common cause of chronic neurological disability in young adults^1,2. Here, to provide insight into the potential mechanisms involved in progression, we conducted a genome-wide association study of the age-related MS severity score in 12,584 cases and replicated our findings in a further 9,805 cases. We identified a significant association with rs10191329 in the DYSF-ZNF638 locus, the risk allele of which is associated with a shortening in the median time to requiring a walking aid of a median of 3.7 years in homozygous carriers and with increased brainstem and cortical pathology in brain tissue. We also identified suggestive association with rs149097173 in the DNM3-PIGC locus and significant heritability enrichment in CNS tissues. Mendelian randomization analyses suggested a potential protective role for higher educational attainment. In contrast to immune-driven susceptibility³, these findings suggest a key role for CNS resilience and potentially neurocognitive reserve in determining outcome in MS.

Extended Data Fig. 1 |. Demographic characteristics by population and center.

a, Discovery population (n = 12,584). b, Replication population (n = 9,805). Bars represent the proportion of patients in each category. Centers are ordered as in the box plot legend (bottom right subpanel). Box plots show median, first, and third quartiles; whiskers represent the smallest and largest values within 1.5-times the interquartile range; outliers are depicted as dots. The countries corresponding to the abbreviations in the box plot legend are shown in Supplementary Table 1. ARMSS, age-related multiple sclerosis severity; EDSS, expanded disability status scale; Primary prog., primary progressive; yrs, years.

Extended Data Fig. 2 |. Principal component analysis of the discovery and replication populations.

MS cases were recruited from 13 countries for the discovery (a) and 8 for the replication (b). After removing population outliers, all remaining cases were of European ancestry. The first two principal components respectively captured the north-to-south and east-to-west gradients of European genetic structure. US and Canadian participants overlapped with those from other countries. Based on self-reported ancestry, East European and Ashkenazi Jewish individuals constituted the majority of the predominantly US subcluster located at the bottom right of the discovery population (a). The scree plots for our principal component analysis in the discovery (c) and replication (d) populations confirm that the first few principal components capture most of the variance attributable to the minimal population structure remaining after quality control.

Extended Data Fig. 3 |. Replication of MS severity variants by center.

a, Genome-wide significant lead variant rs10191329. b, Suggestive lead variant rs149097173. Forest plots show successful replication of the two variants with minimal heterogeneity between centers as indicated by the Cochran’s Q and I² statistics (n = 9,805 participants). ARMSS scores are rank-based inverse-normal transformed. Error bars represent 95% CIs. ARMSS, age-related multiple sclerosis severity; CI, confidence interval.

Extended Data Fig. 4 |. Association of rs149097173 with longitudinal disability outcomes.

a, Adjusted mean EDSS scores over time by carrier status for rs149097173 predicted from LMM analysis. Shaded ribbons indicate the standard error of the mean over time; P value from LMM. b, Covariate-adjusted cumulative incidence of 24-week confirmed disability worsening for the same groups of individuals. c, Covariate-adjusted cumulative incidence of requiring a walking aid; carriers had a 2.2-year shorter median time to require a walking aid. HR and two-sided P values were obtained from Cox proportional hazards models using imputed allele dosage (b–c; Methods). Results were not significant after adjusting for multiple testing across two variants (see Fig. 3 for rs10191329 associations) and three outcomes (P < 0.05/6), although the latter are not expected to be independent. CI, confidence interval; HR, hazard ratio.

Extended Data Fig. 5 |. Tissue expression for nominated MS severity genes.

Gene expression profiles were obtained from GTEx (version 8). Transcripts were collapsed to the gene level and expressed in natural log-transformed transcript per million (TPM) units. DYSF, ZNF638, DNM3 and PIGC are expressed in the brain. Box plots show median, first, and third quartiles; whiskers represent the smallest and largest values within 1.5-times the interquartile range; outliers are depicted as dots. Bold x-axis labels identify CNS tissues. Colors represent tissue types as defined in GTEx.

Extended Data Fig. 6 |. Cell type expression profiles for nominated MS severity genes.

Single-cell RNA sequencing data from 25 human tissues and peripheral blood mononuclear cells were obtained from the Human Protein Atlas. Transcript expression levels were summarized per gene and reported as average normalized transcripts per million (nTPM) in 76 cell types. Asterisks mark cell type specificity for the gene, defined as at least fourfold higher expression in a cell type compared to the mean of others. We note that three of the genes show specificity for oligodendrocyte lineage cells. PIGC expression in brain neuronal and glial cells, missing here, is demonstrated in Extended Data Fig. 8. Colors represent cell type categories; bold x-axis labels identify neuronal and glial cell categories.

Extended Data Fig. 7 |. Cell type expression for PIGC in brain white matter tissue.

Single nuclear RNA expression from 4 progressive MS patients and 5 non-neurological controls confirms PIGC expression in neuronal and glial cells including oligodendrocyte lineage cells. COPs, committed oligodendrocyte precursors; ImOLGs, immune oligodendroglia; Oligo, oligodendrocyte; OPCs, oligodendrocyte precursor cells; Vasc, vascular.

Extended Data Fig. 8 |. Genetic correlations with MS severity.

Shared genetic contribution obtained from cross-trait LDSC. Colors correspond to genetic correlation (r_g) estimates (blue, negative; red, positive). An asterisk indicates a correlation that is significantly different from zero, based on two-sided P values calculated using LDSC (*FDR < 0.05, **FDR < 0.01). Full results are in Supplementary Table 17. Aging-GIP1 was constructed using principal component analysis to capture GWASs of healthspan, father lifespan, mother lifespan, longevity, frailty, and self-rated health.

Extended Data Fig. 9 |. Association of individual MS susceptibility variants (n = 209) with longitudinal disability outcomes.

a, Distribution of P values from adjusted LMM analysis of EDSS change across all study visits. Distribution of two-sided P values from adjusted Cox proportional hazards analyses of (b) time to 24-week confirmed disability worsening and (c) time to require a walking aid. The dashed orange line represents the Bonferroni-corrected significance threshold adjusted for the number of susceptibility variants. d, Venn diagram of nominal associations (P_unadjusted < 0.05) between individual MS susceptibility variants and all disability outcomes considered; no variant showed consistent association across three or more outcomes. The labels in this panel correspond to the following outcomes: ARMSS, association with ARMSS scores following rank-based inverse normal transformation; Disability worsening, time to 24-week confirmed disability worsening; Walking aid, time to require a walking aid (EDSS 6.0); EDSS rate, rate of EDSS change across all study visits.

Extended Data Fig. 10 |. MS susceptibility PGS and longitudinal disability outcomes.

a, Adjusted mean EDSS scores over time by PGS quartile predicted from LMM analysis. Shaded ribbons indicate the standard error of the mean over time; P value from LMM. b, Covariate-adjusted cumulative incidence of 24-week confirmed disability worsening comparing individuals in the highest versus those in the lowest quartile of MS susceptibility PGS. c, Covariate-adjusted cumulative incidence of requiring a walking aid for the same groups of individuals. HR and two-sided P values were obtained from Cox proportional hazards models using imputed allele dosage (b–c; Methods). Across all analyses, the MS susceptibility PGS had no influence on longitudinal outcomes.

Fig. 1 |. Tissue and cell type heritability enrichment.

a, MS susceptibility from previous meta-analysis. b, MS severity from this study. While susceptibility associations display strong immunological lymphoid and myeloid enrichment, our analysis of MS severity uncovered significant enrichment exclusively in CNS tissues. Each point represents one of 205 tissues and cell types, grouped by color into 9 categories. Large circles are significant at a false discovery rate cutoff of 0.05 (dotted line). Full results including tissue and cell type labels are provided in Supplementary Tables 7 and 8.

Fig. 2 |. Within-cases GWAS identifies a novel locus associated with MS severity.

a, Genome-wide association statistics obtained by linear regression of ARMSS scores. The −log₁₀(P) are plotted against chromosomal position. The horizontal dashed line corresponds to the genome-wide significant threshold (P < 5×10⁻⁸) and the horizontal dotted line reflects the threshold for suggestive association (P < 5×10⁻⁶). The bold label indicates the lead genome-wide significant and replicated variant. Variants labeled in gray were not replicated. b, Locus Zoom plot for rs10191329 (DYSF-ZNF638 locus). c, Locus Zoom plot for rs149097173 (DNM3-PIGC locus). Top, −log₁₀(P) from the GWAS for variants at each locus (left y-axis) with the recombination rate indicated by the blue line (right y-axis); bottom, gene positions. Colors represent LD (r² values) with the lead variant. LD, linkage disequilibrium.

Fig. 3 |. MS severity variant accelerates disability accumulation in longitudinal analysis.

a, Adjusted mean EDSS scores over time predicted from LMM analysis showed faster disability worsening in rs10191329 risk allele carriers. Shaded ribbons indicate the standard error of the mean over time; P value from LMM. b, Covariate-adjusted cumulative incidence of 24-week confirmed disability worsening in MS patients based on rs10191329 genotype. Similar to MS clinical trials, worsening was defined as an increase in EDSS by 1.0 if the baseline score was < 5.5 and by 0.5 if the baseline was ≥ 5.5. c, Covariate-adjusted cumulative incidence of requiring a walking aid for the same lead variant. Homozygous carriers had a 3.7-year shorter median time to require a walking aid. HR and two-sided P values were obtained from Cox proportional hazards models using imputed allele dosage (b–c; Methods). Left-censoring of participants with EDSS ≥ 6.0 at study entry resulted in different sample sizes for genotype groups in the time to walking aid analysis. CI, confidence interval; HR, hazard ratio.

Fig. 4 |. Cortical lesion rate and brainstem lesion count are elevated in homozygous rs10191329 risk allele carriers.

a, Schematic representation of tissue sampling locations (created with BioRender.com). Demyelinating lesions were quantified on a brainstem section dissected in a consistent manner across individuals. Cortical lesions were identified on supratentorial tissue blocks targeted to macroscopic or MRI-visible MS lesions. b, Brain tissue section immunostained for the proteolipid protein marker of myelin (brown). A subpial cortical lesion characterized by loss of myelin is marked by an asterisk and delineated by the dotted white line. The solid white line separates normal-appearing gray matter (sparse brown) from white matter (dense brown). c, A lesion spanning gray and white matter in the brainstem of the same donor, marked by an asterisk and delineated from normal-appearing tissue by the dotted white line. Black scale bars indicate 0.5 mm. The donor was an A allele homozygote for rs10191329. d, Box plots show median, first, and third quartiles; whiskers represent the smallest and largest values within 1.5-times the interquartile range; outliers are depicted as dots. Two-sided P values were obtained from generalized linear models comparing lesion count in the cortex (offset by the relevant number of tissue blocks; n = 174 donors) and brainstem (n = 181 donors) across genotype groups adjusting for covariates; significant differences are marked with an asterisk. The displayed cortical lesion rate was calculated by dividing the number of lesions by the number of tissue blocks containing cortex.

Fig. 5 |. Association of MS severity with educational attainment and smoking.

a, MR estimates for the effect of years of education (n = 765,283), lifetime smoking index (n = 462,690), body mass index (n = 681,275), and 25-hydroxyvitamin D (n = 441,291) on ARMSS scores; the lighter color represents nonsignificant results. b, Similarly, adjusted polygenic risk score (n = 12,584) and observational analyses of two MS cohorts (n = 2,878 and 5,228) demonstrated reduced MS severity with higher years of education in linear regression models. This effect persisted following adjustment for smoking and income. c, Mean ARMSS scores decreased with higher PGS_EDU percentile. d, Similarly, higher percentile of recorded years of education associated with lower mean ARMSS scores in the EIMS cohort. e, Mean ARMSS scores decreased with higher percentile years of education in the GEMS cohort. P values were obtained from a regression of ARMSS scores on PGS_EDU (c) or years of education (d-e), adjusted for baseline covariates. In the MR and observational analyses, point estimates (squares) reflect a 1-year increase in education, while polygenic score estimates are per standard deviation score increase. ARMSS, age-related multiple sclerosis severity score (rank-based inverse-normal transformed); IVW, inverse-variance weighted; PGS_EDU, polygenic score for years of education; RAPS, robust adjusted profile score; WM, weighted median.