Genetic risk and a primary role for cell-mediated immune mechanisms in multiple sclerosis.

Multiple sclerosis is a common disease of the central nervous system in which the interplay between inflammatory and neurodegenerative processes typically results in intermittent neurological disturbance followed by progressive accumulation of disability. Epidemiological studies have shown that genetic factors are primarily responsible for the substantially increased frequency of the disease seen in the relatives of affected individuals, and systematic attempts to identify linkage in multiplex families have confirmed that variation within the major histocompatibility complex (MHC) exerts the greatest individual effect on risk. Modestly powered genome-wide association studies (GWAS) have enabled more than 20 additional risk loci to be identified and have shown that multiple variants exerting modest individual effects have a key role in disease susceptibility. Most of the genetic architecture underlying susceptibility to the disease remains to be defined and is anticipated to require the analysis of sample sizes that are beyond the numbers currently available to individual research groups. In a collaborative GWAS involving 9,772 cases of European descent collected by 23 research groups working in 15 different countries, we have replicated almost all of the previously suggested associations and identified at least a further 29 novel susceptibility loci. Within the MHC we have refined the identity of the HLA-DRB1 risk alleles and confirmed that variation in the HLA-A gene underlies the independent protective effect attributable to the class I region. Immunologically relevant genes are significantly overrepresented among those mapping close to the identified loci and particularly implicate T-helper-cell differentiation in the pathogenesis of multiple sclerosis.

We performed a large GWAS as part of the Wellcome Trust Case Control Consortium 2 (WTCCC2) project. Cases were recruited through the International Multiple Sclerosis Genetics Consortium (IMSGC) and compared with the WTCCC2 common control set[12,13] supplemented by data from the control arms of existing GWAS. We introduced a number of novel quality control (QC) methods for processing these datasets (see Supplementary Information), which ultimately provided reliable information from 9772 cases and 17376 controls (see Figure 1A). Following single nucleotide polymorphism (SNP) based QC, data from 441547 autosomal SNPs, common to all internally and externally generated datasets, were available for analysis.

Figure 1

Distribution of cases and controls. All cases and controls were drawn from populations with European ancestry; cases from 15 countries and controls from 8. A: numbers of case (red) and control (black) samples from each country. B: The projection of samples onto the first two principal components of genetic variation, with cases shown on the left and controls on the right. The axes are orientated to approximate to the geography, and samples are colour coded as indicated in the legend. We genotyped the cases (9772) and some Swedish controls (527) using the Illumina Human 660-Quad platform, and the UK controls (5175, the WTCCC2 common control set[12,13]) using the Illumina 1.2M platform. All other controls were genotyped externally using various Illumina genotyping systems (see Supplementary Information).

The multi-population nature of our study (Figure 1 A and B) afforded an opportunity to assess various published approaches for controlling the potential confounding effects of population structure, several of which (in the event) proved unhelpful (see Supplementary Information). Whilst not common in primary GWAS undertaken to date, the challenge of combining data across populations, in contexts where not all case samples have controls available from the same population (thus precluding standard meta-analytical techniques), may become more routine as study sizes increase.

We attempted analyses of the non-United Kingdom (UK) data with the now widespread technique of using principal components as covariates to correct for structure. However, even use of all seven top principal components which captured genome-wide effects in our data resulted in an unacceptably high genomic inflation: for example, the genomic control factor (λ)[14] was λ = 1.2. We tried to reduce the genomic inflation by discarding the case samples that seemed least well matched to control sets. Removal of half the available cases in this fashion only reduced λ to 1.1. In another approach to handling structure, statistical clustering algorithms were successful in identifying subgroups of the data within which cases and controls appeared well-matched for ancestry (See Supplementary figure S17). However tests within these subgroups combined via fixed-effects meta-analysis also yielded unacceptably high genomic inflation (λ >1.4) in an analysis with seven matched sub-groups of cases and controls. Finally, we applied a novel variance components method (similar to Kang et al.[15]), separately to the UK and non-UK datasets, that explicitly accounts for correlations among the phenotypes of individuals resulting from relatedness, allowing us to deal successfully with all sources of structure in our samples (see Supplementary Information for details of the linear mixed model we used). For example, the genomic inflation was reduced to λ = 0.995 in the UK and 1.016 in the non-UK data (see also Supplementary Information). After fixed effects meta-analysis of the results from the UK and non-UK datasets, the inflation factor was λ = 1.045. We adopted this approach for all subsequent non-MHC association analyses.

Outside the MHC we identified 95 distinct regions having at least one SNP associated with multiple sclerosis at pGWAS <1×10−4.5; in six of these 95 regions conditional analysis revealed an additional SNP showing association to the same locus (one locus containing two such SNPs). In total we took all 102 SNPs forward to replication, which we performed using data from previously reported multiple sclerosis GWAS[8,9] and the iControl database (excluding any WTCCC controls previously used in these studies). In total, the replication analysis included data from 4218 cases and 7296 controls. These were considered in six independent strata after which results were combined through a fixed effects meta-analysis. For 98 of the 102 SNPs, the same allele was over-represented in cases compared to controls. Twenty three of the 26 previously known or strongly suggested multiple sclerosis associated loci were replicated in our primary GWAS with pGWAS <1×10−3. Our GWAS and replication also revealed another 29 novel associated regions (defined as having pGWAS <1×10−4.5, one-sided pReplication <0.05, and pCombined <5×10−8), and a further 5 regions with strong evidence for association (with pGWAS <1×10−4.5, one sided pReplication <0.05, and pCombined <5×10−7). In one previously reported locus and two novel loci, additional SNPs were identified as being conditionally important in explaining risk. Just over one third of the identified loci overlap with regions already confirmed as associated with at least one other autoimmune disease (according to the GWAS catalog, http://www.genome.gov/gwastudies/). Results both for the previously established and novel loci are shown in Figure 2 and Supplementary Tables 1-3; and details of all 102 SNPs taken to replication are available in the Supplementary Data file.

Figure 2

Regions of the genome showing association to multiple sclerosis. Columns from left to right: evidence for association from the linear mixed model analysis of the discovery data (thresholded at −log10(p) = 12). Non-MHC regions containing associated SNPs are shown in red and are labelled with the rsID (bold for newly identified loci, black for strong evidence, grey for previously reported) and risk allele of the most significant SNP. * indicates that the locus contains a secondary SNP signal. Odds ratio and 95% confidence intervals estimated from the meta-analysis of the discovery and replication data (+ indicates estimates for previously-known loci from discovery data only). Risk allele frequency estimates in each of the control populations used in the study (each is shown as a vertical bar on a scale from 0 to 1 going left to right). For each region of association the number of genes is reported, and where non-zero a candidate gene is given. Black dots indicate that the candidate gene is physically the nearest gene included in the “immune system process” GO term. When the most-significant SNP tags a SNP predicted to have an impact on the function of the candidate gene this is indicated. Where such a SNP exists, the gene involved is selected as the candidate gene; otherwise the nearest gene is selected unless there are strong biological reasons for a different choice. The final column indicates SNPs which are correlated (r > 0.1) with SNPs reported to be associated with other autoimmune (AI) diseases (abbreviations: RA = Rheumatoid arthritis; CeD = Celiac disease; UC = Ulcerative colitis; CrD = Crohns’s disease; T1D = Type 1 diabetes; PS = Psoriasis). An interactive version of the figure is available at www.well.ox.ac.uk/wtccc2/ms.

In order to assess objectively the collective evidence across the associated regions for particular classes of genes, we performed statistical analyses to look for enrichment of genes with similar function. We first identified the nearest gene to the lead SNP in each of the (52) regions of association and used the Gene Ontology (GO) database[16] to define sets of functionally related genes (GO terms). We then tested whether the set of nearest-genes was enriched for particular GO terms using Fisher’s exact test. The GO terms having the most significant enrichment include genes linked to lymphocyte function (p =3.2×10−11, OR = 35.96) and in particular those with a role in T cell activation and proliferation (p = 1.85×10−9, OR = 40.85). These are representative of a larger group associated with various components of the GO ‘immune system process’ (p = 8.6×10−11, OR = 9.12). A similar analysis based on all genes in or near association regions showed similar enrichment, as did independent analyses based on nearest-gene or all genes in our next tier of signals, the 42 regions taken to replication but not meeting the thresholds above for association (see Supplementary file.) Although GO immune system genes only account for 7% of human genes, in 30% of our association regions the nearest gene to the lead SNP is an immune system gene. As an illustration, Figure 3 shows a schematic of genes involved in the T helper cell differentiation pathway; a striking number show strong evidence for association with multiple sclerosis particularly those acting as cell surface receptors. We infer from this pathway analysis of our GWAS signals that specific classes of immune system genes are especially important in the pathogenesis of multiple sclerosis.

Figure 3

Graphical representation of the T helper cell differentiation pathway. The figure is derived from an image generated by Ingenuity Pathway Analysis (IPA) software version 8.8 (Ingenuity Systems, Inc., Redwood City, CA, USA). Alpha-numeric labels indicate the individual genes and gene complexes (nodes) included in the pathway (note some are included more than once). Coloured nodes are those containing a gene implicated by proximity to a SNP showing evidence of association. Red: in bold or grey in Figure 2 (plus MHC class II region and TNFα); Orange: other loci in Figure 2 or discovery P value < 1×10−4.5 and consistent replication data. Yellow: Discovery P value < 1×10−3. Other molecules (proteins, vitamins etc.) may also be of relevance in these processes but are not included here as they are not currently listed as being part of this particular pathway in the IPA database.

Our screen not only implicates a multitude of genes coding for cytokine pathway (CXCR5, IL2RA, IL7R, IL7, IL12RB1, IL22RA2, IL12A, IL12B, IRF8, TNFRSF1A, TNFRSF14, TNFSF14), co-stimulatory (CD37, CD40, CD58, CD80, CD86, CLECL1) and signal transduction (CBLB, GPR65, MALT1, RGS1, STAT3, TAGAP, TYK2) molecules of immunological relevance, but also relates to previously reported environmental risk factors such as Vitamin D[9,17] (CYP27B1, CYP24A1) and therapies for multiple sclerosis including Natalizumab[18] (VCAM1) and Daclizumab[19] (IL2RA). There is a relative absence of genes relevant to potential pathways for neurodegeneration independent of inflammation (GALC, KIF21B).

To refine our understanding of the MHC associations in multiple sclerosis we imputed classical Human Leukocyte Antigen (HLA) types at six loci (A, B, C, DQA1, DQB1 and DRB1)[20] and analysed these alongside the SNPs (see Supplementary Information for validation; at alleles responsible for the major signals described below, estimated specificity was at least 0.99 and sensitivity was at least 0.98, except for DRB1*13:03, where it was 0.88). Primary discovery was focused on the UK cohort with candidate signals being validated through support from additional case-control cohorts. Because of the extensive linkage disequilibrium within the MHC, we identified associated alleles in a stepwise manner, selecting the most strongly associated to include in a general model, in turn, if pUK <10−4 and pcombined <10−9 (Supplementary Information). At each stage we explored possible interactions and departures from the simple model in which risk increases multiplicatively with each additional copy of the relevant allele (additive increase on the log-odds scale) within the logistic risk framework.

By this approach we found that DRB1*15:01 has the strongest association with multiple sclerosis amongst all classical and SNP alleles, with a consistent effect between cohorts (p <1×10−320: Figure 4A). The data are consistent with an additive effect on the log-odds scale for each additional allele. Conditioning on DRB1*15:01, we confirmed the presence of a protective Class I allele and identified the signal as being driven by HLA-A*02:01 (as previously suggested),[21] with a consistent effect size across cohorts (p = 9.1×10−23: Figure 4A). Again, we found no strong evidence for departure from additivity on the log-odds scale or statistical interaction with DRB1*15:01. Conditioning on both DRB1*15:01 and A*02:01 revealed additional risk associated with the strongly linked alleles DRB1*03:01 and DQB1*02:01 (p = 3.6×10−10: Figure 4A; note that we cannot separate these alleles but for simplicity refer only to DRB1*03:01 below). Further conditioning identified an additional DRB1 risk allele DRB1*13:03 (p = 1.3×10−11: Figure 4A). Although no other classical alleles meet the above criteria, we did observe several SNPs providing independent signals, the strongest coming from rs9277535_G (combined OR 1.28, p = 2.2×10−22), an allele known to be in linkage disequilibrium with DPB1*03:01 (r2 = 0.37).[22]

Figure 4

Results for the main MHC alleles. A: Forest plots for each of the primary HLA alleles (HLA-A*02:01, DRB1*15:01, DRB1*03:01 and DRB1*13:03) showing consistency of effect across the populations and combined OR of 0.73, 3.1, 1.26 and 2.4 respectively (whiskers indicate 95% confidence intervals). B: The genealogical trees estimated for DRB1 (top) and HLA-A (bottom). These trees were constructed using classical HLA and SNP typing data available from the HapMap CEU haplotype data. Each left hand branch of the tree terminates on a set of haplotypes carrying a particular HLA allele. The coloured dots indicate the mostly likely locations for a disease-associated mutation as predicted by the GENECLUSTER program[23]. In the DRB1 tree, the blue dot captures a risk effect attributable to all haplotypes carrying the *15:01 allele. The green dot captures a risk effect carried by all haplotypes carrying the *03:01 allele and the red dot captures a risk effect on haplotypes carrying *13:03 or *08:01. In the HLA-A plot, the orange dot is a protective mutation lying at the root of all *02:01, *02:05, *02:06 and *68:01 alleles. The blue dot in brackets denotes a branch containing those *03:01 haplotypes that also carry DRB1*15:01; the GENECLUSTER prediction here is thus a reflection, due to linkage disequilibrium of the risk attributable to DRB1*15:01. The terminal branches are labelled with the allele carried by the haplotype and its frequency.

Analysis of the MHC SNP data using a genealogical method (GENECLUSTER)[23] offers an alternate means of relating our results to classical HLA alleles that provides additional insight into the underlying genetic architecture (see Supplementary Information). Figure 4B shows genealogical trees relating the classical alleles at DRB1 and HLA-A, together with the estimated evolutionary position of the mutations predicted by GENECLUSTER as most completely modelling the association. At HLA-DRB1, three mutations are predicted, each of which implicates a clade of haplotypes carrying particular DRB1 alleles. All of the DRB1 alleles we have shown to be independently associated are included in these clades, each corresponding to a particular mutation. In addition, the analysis also explains why those haplotypes carrying the *08:01 allele have previously been shown to increase risk[24,25] since they carry the same mutation as those bearing *13:03. At HLA-A, the predicted protective mutation is also concordant with our regression analysis of classical alleles in implicating *02:01 but, in addition, predicts that *68:01, *02:05, and *02:06 carry the same protective allele. All of these secondary predictions (increased risk from DRB1*08:01 and protection from HLA-A*68:01, *02:05, and *02:06) are supported in our regression analysis of classical alleles but the power to detect them in the primary analyses is limited because each allele occurs at a very low frequency.

We found no evidence for genetic associations with clinical course, severity of disease or month of birth, and no evidence of interaction with gender or DRB1*15:01 in any part of the genome (see Supplementary Information). However, analysis with respect to age at onset replicated the previously suggested association with the DRB1*15:01 allele.[26] Although no other part of the genome contained individual SNPs showing strong evidence for association, risk alleles determining susceptibility are collectively more closely associated with age at onset than expected by chance, suggesting that individual genetic susceptibility is inversely correlated with age at onset.

Our GWAS - large for any complex trait having a prevalence of 1:1000 and involving diverse populations of European descent - has identified 29 novel susceptibility loci. Four mutations, one from Class I and three from Class II, with effects modelled in a simple multiplicative manner within and across loci are sufficient to account for most of the risk attributable to the MHC (see Supplementary file). Although our data do not address the issue of which components within the nervous system are initially damaged by the inflammatory response the over-representation of genes that influence T cell maturation provides independent and compelling evidence that the critical disease mechanisms primarily involve immune dysregulation.

More generally, our study reinforces the view that the GWAS design, combined with very large experimental sample sizes and careful statistical analysis, provides valuable insights into the genetic architecture of common complex diseases. Here, this approach has identified many associated genetic variants close to genes, which are both individually interesting and collectively illuminate the roles of key biological pathways. It also provides indirect evidence that many more common variants of small effect contribute to genetic susceptibility for multiple sclerosis. Simple models, in which the previously-known and newly-identified variants affect risk multiplicatively, both within and across loci, explain a meaningful proportion (~20%, see Supplementary Information) of genetic risk for the disease. Important challenges lie ahead, in understanding overlap between the genetic basis for susceptibility in the context of different autoimmune diseases, and in uncovering the functional mechanisms underlying these associations.

Materials and Methods

Details of case ascertainment, processing and genotyping, together with sample and genotyping quality control are provided in the Supplementary Information. Statistical methods developed for testing the reliability of externally generated data sets, detecting samples with non-European ancestry, correcting for structure, classical HLA imputation and meta-analysis are also outlined in the Supplementary Information. Results for all scans and all reported loci are described in detail in the Supplementary Information.