A knowledge-driven interaction analysis reveals potential neurodegenerative mechanism of multiple sclerosis susceptibility.

Gene-gene interactions are proposed as an important component of the genetic architecture of complex diseases, and are just beginning to be evaluated in the context of genome-wide association studies (GWAS). In addition to detecting epistasis, a benefit to interaction analysis is that it also increases power to detect weak main effects. We conducted a knowledge-driven interaction analysis of a GWAS of 931 multiple sclerosis (MS) trios to discover gene-gene interactions within established biological contexts. We identify heterogeneous signals, including a gene-gene interaction between CHRM3 (muscarinic cholinergic receptor 3) and MYLK (myosin light-chain kinase) (joint P=0.0002), an interaction between two phospholipase C-β isoforms, PLCβ1 and PLCβ4 (joint P=0.0098), and a modest interaction between ACTN1 (actinin alpha 1) and MYH9 (myosin heavy chain 9) (joint P=0.0326), all localized to calcium-signaled cytoskeletal regulation. Furthermore, we discover a main effect (joint P=5.2E-5) previously unidentified by single-locus analysis within another related gene, SCIN (scinderin), a calcium-binding cytoskeleton regulatory protein. This work illustrates that knowledge-driven interaction analysis of GWAS data is a feasible approach to identify new genetic effects. The results of this study are among the first gene-gene interactions and non-immune susceptibility loci for MS. Further, the implicated genes cluster within inter-related biological mechanisms that suggest a neurodegenerative component to MS.

Introduction

Multiple sclerosis (MS) is a complex autoimmune disorder characterized by demyelination and scar tissue formation within the central nervous system. Axonal loss and progressive neurodegeneration leads to impaired neurological function and diminished quality of life. The MHC region of chromosome 6 is consistently associated with MS risk1, 2, and several additional immunological and inflammatory loci have recently been implicated in MS susceptibility. However, all known susceptibility loci combined account for far less than 50% of the estimated heritability.

Epistasis or gene-gene interaction has been promoted as an important part of complex disease etiology, due to the monumental complexity of biological systems. In addition, when explicitly modeled, it can increase power to detect the independent effects of susceptibility loci 3, 4. However, a key challenge of interaction analysis is the biological interpretation of statistical results and the task of resolving functional relationships between variants from multiple loci across the genome. With these challenges in mind, we assessed gene-gene interactions within established biological contexts based on multiple knowledge sources, such as pathway and ontology databases. This analysis implicates several new functionally related MS risk loci, including three interaction effects centered on calcium signaling, representing a neurodegenerative genetic component to MS etiology.

Results and Discussion

For our analysis, we examined gene-gene interactions in a screening dataset followed by three validation sets. The screening phase was conducted in two stages to maximize information gain from the analysis. Using the trio study design of our screen dataset, we examined the transmission of alleles to affected offspring. This stage revealed 5,463 models (consisting of 5,965 SNPs) with MF p < 0.001 and LR p < 0.001, indicating an effect from simultaneous over-transmission of alleles at two separate loci. The majority of these models contain distinct SNPs, however some are recurrent across multiple models. The top recurring SNP was involved in 35 models, and most occurred in less than 5 models. Also notably the majority of models do not contain the known main effect of the MHC region, with only 474 of 5,463 containing markers from chromosome 6 where the MHC resides. The joint transmission of alleles to offspring could also be due to chromosomal linkage if the two alleles are in close enough proximity that recombination events between the two loci are infrequent. To minimize false-positives due to this possibility and to provide additional evidence for interaction, we compared the probands from the 931 families with 2,950 controls from the Wellcome Trust Case Control Consortium (WTCCC)5. Using standard logistic regression on the proband/control set, we reduced the 5,463 models to 326 models (MF p < 0.001 and LR p < 0.001) containing 469 SNPs. For all validation sets, we established α = 0.05 for both MF and LR tests as our replication criteria.

In validation set I, 306 of the 326 models contain all SNPs available after QC procedures and were evaluated. Twenty multi-locus models had significant model fit (MF p < 0.05) and direction of coefficients consistent with the screening dataset, representing nine gene-gene combinations (supplemental table 1).

While this analysis approach incorporated prior biological knowledge, it was somewhat unexpected that of hundreds of thousands of gene-gene relationships considered, four of the nine models with statistically consistent results point to two interconnected biological mechanisms -- a calcium-signaled change in cytoskeleton dynamics. Of these four models, two had a significant interaction component (LR p < 0.05). To examine and validate the role of this biological mechanism in MS susceptibility, we explored the four functionally related models (table 1) in two additional datasets, validation set II and III. Validation set III was genotyped with a different platform from the other studies, so surrogate SNPs were selected with close physical proximity and high r-squared based on the Hapmap CEU data.

Table 1

Significant Models from Screen and Validation Set I localized to calcium signaling and cytoskeleton regulation

#	Locus 1			Locus 2			Screen Trio Conditional LR		Screen Proband/Control LR		Validation Set I
	Chr	Gene	SNP	Chr	Gene	SNP	Model Fit	Interaction	Model Fit	Interaction	Model Fit	Interaction
1	7	SCIN	rs2240571	15	CYFIP1	rs8025779	3.75E-04	1.51E-04	0.0001	0.0001	0.0049	0.3565
2	14	ACTN1	rs17106421	22	MYH9	rs1009150	8.93E-04	6.38E-05	0.0001	0.0001	0.0082	0.0952
3	1	CHRM3	rs528011	3	MYLK	rs4677905	5.57E-04	3.74E-05	0.0005	0.0001	0.0235	0.0025
4	20	PLCB4	rs4816129	20	PLCB1	rs6516415	9.23E-04	8.50E-05	0.0008	0.0009	0.0443	0.0095

The results from validation sets II and III indicate heterogeneous weak effects from variants within this biological mechanism (table 2). In validation set II, model 1 had a non-significant model fit (MF = 0.0595), and there was no evidence of interaction (LR = 0.7664). In validation set III, model 1 also had a non-significant model fit (MF = 0.0748) but with a significant interaction term (LR = 0.0366). Models 2, 3, and 4 had non-significant p-values (MF and LR) in both validation sets, and while model 1 was not significant, it did show consistently low p-values (p < 0.1) across all evaluated datasets.

Table 2

Details for models related to calcium signaling and cytoskeleton regulation evaluated in screen and validation sets.

						Interaction Model					Main Effect Models
	#	SNP1	MAF	SNP2	MAF	MF p	LR p	Coeff 1	Coeff 2	Coeff 1X2	SNP1 p	OR	SNP2 p	OR
screen	1	rs2240571	0.49	rs8025779	0.37	0.0001	0.0001	1.64 (1.3,2.08)	1.55 (1.26,1.91)	0.72 (0.62,0.85)	0.2014	1.07 (0.96,1.19)	0.1092	1.10 (0.98,1.22)
	2	rs17106421	0.05	rs1009150	0.28	0.0001	0.0001	5.25 (2.32,11.89)	6.71 (2.66,16.91)	0.36 (0.23,0.59)	0.3859	1.11 (0.87,1.41)	0.5390	0.96 (0.86,1.08)
	3	rs528011	0.44	rs4677905	0.25	0.0005	0.0001	1.55 (1.16,2.07)	1.53 (1.21,1.94)	0.71 (0.59,0.84)	0.0967	0.92 (0.82,1.02)	0.5363	1.04 (0.92,1.18)
	4	rs4816129	0.35	rs6516415	0.34	0.0008	0.0009	0.73 (0.58,0.92)	0.62 (0.49,0.78)	1.31 (1.12,1.54)	0.8271	1.01 (0.91,1.13)	0.0245	0.88 (0.79,0.98)
validationI	1	rs2240571	0.47	rs8025779	0.36	0.0049	0.3565	1.38 (1.06,1.79)	1.10 (0.87,1.38)	0.92 (0.77,1.1)	0.0004	1.24 (1.1,1.4)	0.9689	1.00 (0.88,1.13)
	2	rs17106421	0.05	rs1009150	0.31	0.0082	0.0952	1.75 (0.89,3.41)	1.64 (0.72,3.72)	0.70 (0.46,1.07)	0.7015	1.06 (0.8,1.39)	0.0022	0.82 (0.72,0.93)
	3	rs528011	0.45	rs4677905	0.24	0.0230	0.0025	1.69 (1.2,2.38)	1.44 (1.09,1.9)	0.73 (0.59,0.9)	0.5333	1.04 (0.92,1.17)	0.9694	1.00 (0.87,1.15)
	4	rs4816129	0.36	rs6516415	0.34	0.0443	0.0090	0.70 (0.53,0.92)	0.79 (0.61,1.02)	1.28 (1.06,1.53)	0.5965	0.97 (0.86,1.09)	0.3274	1.06 (0.94,1.21)
validationII	1	rs2240571	0.48	rs8025779	0.36	0.0595	0.7664	1.09 (0.92,1.30)	1.00 (0.86,1.16)	1.02 (0.90,1.15)	0.0073	1.12 (1.03,1.22)	0.6782	1.02 (0.94,1.11)
	2	rs17106421	0.05	rs1009150	0.30	0.2102	0.5750	0.99 (0.63,1.55)	0.93 (0.53,1.62)	1.09 (0.81,1.46)	0.2884	1.11 (0.92,1.34)	0.0862	1.08 (0.99,1.19)
	3	rs528011	0.46	rs4677905	0.23	0.6175	0.2773	1.10 (0.87,1.39)	1.11 (0.92,1.34)	0.92 (0.8,1.07)	0.4984	0.97 (0.89,1.06)	0.6282	1.02 (0.93,1.13)
	4	rs4816129	0.36	rs6516415	0.35	0.9794	0.6993	0.96 (0.79,1.16)	0.97 (0.8,1.16)	1.03 (0.9,1.17)	0.8517	0.99 (0.91,1.08)	0.9564	1.00 (0.91,1.09)
validationIII	1	rs3735222	0.42	rs7175101	0.49	0.0748	0.0366	1.19 (0.93,1.51)	1.33 (1.06,1.64)	0.81 (0.67,0.99)	0.5689	0.96 (0.83,1.10)	0.1348	1.11 (0.97,1.27)
	2	rs723113	0.48	rs2157256	0.33	0.4737	0.3686	1.23 (0.91,1.68)	1.14 (0.88,1.48)	0.91 (0.74,1.12)	0.2292	1.08 (0.95,1.24)	0.8449	1.01 (0.88,1.17)
	3	rs1594513	0.24	rs9289225	0.24	0.9461	0.5725	1.11 (0.72,1.71)	1.03 (0.84,1.26)	0.93 (0.71,1.20)	0.8359	0.98 (0.84,1.26)	0.9140	0.99 (0.85,1.16)
	4	rs6118558	0.31	rs6118378	0.28	0.1764	0.4969	0.87 (0.62,1.22)	1.05 (0.76,1.47)	1.08 (0.87,1.34)	0.7899	0.98 (0.85,1.13)	0.0391	1.17 (1.01,1.35)
Joint Typed	1	rs2240571	0.48	rs8025779	0.36	0.0000	0.0677	1.24 (1.1,1.4)	1.13 (1.02,1.26)	0.93 (0.85,1.01)	5.18E-05	1.12 (1.06,1.18)	0.1356	1.04 (0.99,1.11)
	2	rs17106421	0.05	rs1009150	0.30	0.0326	0.0129	1.56 (1.13,2.15)	1.57 (1.06,2.32)	0.78 (0.63,0.95)	0.1998	1.09 (0.96,1.23)	0.2687	0.97 (0.91,1.03)
	3	rs528011	0.45	rs4677905	0.24	0.0002	0.0001	1.29 (1.1,1.5)	1.28 (1.13,1.45)	0.83 (0.75,0.91)	0.1387	0.96 (0.91,1.01)	0.2472	1.04 (0.97,1.11)
	4	rs4816129	0.36	rs6516415	0.34	0.0098	0.0013	0.82 (0.72,0.93)	0.82 (0.73,0.93)	1.15 (1.05,1.25)	0.4858	0.98 (0.93,1.04)	0.4448	0.98 (0.92,1.04)

As models 1 and 2 did not show evidence of interaction in validation set I, and as interaction analysis can enhance the detection of main effects4, we hypothesized that significant models identified in the screening phase could consist of weak effects of independent SNPs. Within validation set I, the significant model fit of models 1 and 2 are driven by main effects of rs1009150 within MYH9 (p = 0.0022) and rs2240571 within SCIN (p = 0.0004). In validation set II, rs2240571 shows a main effect (p = 0.0073). In validation set III, rs6118378 of PLCβ1 (surrogate for rs6516415) has a significant main effect (p = 0.0391).

Combining all datasets with directly typed SNPs (excluding validation set III) for a joint analysis, all models have significant MF p-values (α = 0.05), and all but model 1 have significant LR p-values (α = 0.05). In this joint analysis, SCIN SNP rs2240571 appears as a significant main effect (p = 5.2e-5), and its interaction with CYFIP1 SNP rs8025779 is non-significant (LR p = 0.0677).

Models

Model 1

rs2240571 is located approximately 400 BP upstream from the scinderin (SCIN) gene on chromosome 7. In some datasets, SCIN statistically interacts with cytoplasmic fragile X mental retardation 1 (FMR1) interacting protein 1, or CYFIP1.

SCIN is a calcium-dependent actin severing protein5. In addition to binding calcium, SCIN has PIP2 binding sites that likely play a role in its activity. SCIN mediated disassembly of the cortical actin cytoskeleton may regulate translocation of secretory vesicles to facilitate neurotransmitter release6. SNPs in the SCIN gene curiously did not show a significant main effect from the transmission disequilibrium test (p = 0.1026), Cochran-Mantel-Haenszel test (p = 0.3762), or additive logistic regression analysis (p = 0.2014) of the screen data. In the screen dataset, there was approximately 48% power (α = 0.001) to detect an additive effect of 1.213 (the point estimate for the SNP from the overall dataset), and an effect from rs2240571 was only seen when modeled with rs8025779, a SNP in the CYFIP1 gene. Assessing the joint effect of these SNPs allowed the detection of the weak main effect from the SCIN SNP, either due to a true interaction effect that was oversampled in the screen data, or perhaps due to subtle changes in allele frequency that prevented detection of the interaction in the validation sets7. The effect from SCIN was seen in all datasets, appearing as a main effect in two datasets and as an interaction with CYFIP1 in two others. Though the nature of the interaction effect was inconsistent, the direction of the effect from SCIN was consistent in all models, and in the overall set of samples, SCIN appeared as a main effect.

Model 2

rs17106421 is a SNP 6KB upstream of actinin alpha 1 (ACTN1) located on chromosome 14 and rs1009150 is an intronic SNP in myosin heavy chain 9 (MYH9) located on chromosome 22. ACTN1 and MYH9 function in the formation of actin stress fibers and cytoskeletal contraction. ACTN1 plays a key role in phosphoinositide-3-kinase-induced cytoskeletal reorganization 8, and has brain-specific splicing isoforms 9. Variations in the MYH9 gene have been implicated in a wide variety of disorders including HIV-associated nephropathy, hypertension and other kidney diseases10, and hearing loss11.

Model 3

rs528011 is an intronic SNP located on chromosome 1 in the muscarinic cholinergic receptor 3 (CHRM3) gene between exons 4 and 5. This SNP statistically interacts with rs4677905, an intronic SNP in the myosin-light-chain kinase (MYLK) gene located on chromosome 3. These genes are related by the calcium signaling pathway and also play a role in regulation of the actin cytoskeleton (KEGG pathways: ko04810 & ko04020).

CHRM3 is a G-protein coupled receptor that binds acetylcholine, freeing the G-protein complex to activate phospholipase C (PLC) isoforms, generating inositol 1,45-tripohsphate (IP3). IP3 binds to the IP3-receptor to release Ca2+ ions from intracellular stores. MYLK is activated by this downstream intra-cellular calcium release12. CHRM3 is lost in astrogliotic MS lesions13, and acetylcholinesterase inhibitors reduce the clinical severity of experimental autoimmune encephalomyelitis (EAE), the mouse model of MS14. MYLK mediates myosin II motor activity responsible for actin cytoskelelon contraction, is upregulated during axon regeneration15, and plays a role in axon retraction and regeneration mechanisms16. Differential expression of MYLK is seen in astrocytes from glaucoma patients, and appears to be part of a collection of genes involved in neurodegenerative processes17. As the IL-2 receptor (IL2R) has been implicated in MS by multiple studies18-20, it is noteworthy that intracellular calcium release also triggers multiple downstream events including IL-2 production21.

Model 4

rs4816129 in the intronic region of the phospholipase C beta 4 (PLCβ4) gene is located on chromosome 20 and statistically interacts with intronic SNP rs6516415 in the phospholipase C beta 1 (PLCβ1) gene, located approximately 200 KB upstream. Despite being in close physical proximity, these two SNPs are not in linkage disequilibrium in any of the datasets. PLCβ4 and PLCβ1 function together in multiple KEGG pathways, including the calcium signaling pathway, Wnt signaling, and inositol phosphate metabolism (KEGG: ko04020, ko04310, & ko00562).

PLCβ1 and PLCβ4 are two isozymes in the larger phospholipase-C family22, which hydrolyze phosphatidylinositol 4,5-bisphosphate (PIP2) to produce IP3 and diacylglycerol (DAG). Notably, DAG activates various protein kinase C (PRC) isoforms. Further, PLCβ1 and PLCβ4 are expressed in the central nervous system23, 24, and model systems illustrate a role for both isoforms in proper conduction of nerve signals25.

The genes identified fall among the first susceptibility loci for MS ostensibly involved with the central nervous system and neuron function. Previous studies have identified numerous genes implicated primarily in the autoimmune inflammatory process 26-29. The recent analysis of Baranzini et al. identified general patterns of significance in axon guidance and neurogenesis pathways using the entire GWAS data from Validation set III30, but did not specifically identify any of the genes found in this study. From these results we identify calcium-signaled cytoskeleton regulation as potential neurodegenerative mechanism for MS risk.

Using basic logistic regression procedures, we identified three models that have a non-additive interaction component while significantly contributing to MS risk, and further identify a main effect that was undetected by single-locus analysis of the GWAS data, confirming the principle that interaction analysis can improve detection of weak main effects4. Methodologically, this study is among the first to apply knowledge-based interaction analysis to genome-wide association data30-35, and illustrates that restricting evaluation of two-locus interaction models to those with established biological contexts is a viable strategy.

Materials and Methods

The International Multiple Sclerosis Genetics Consortium (IMSGC) genotyped 931 parent-affected child trios for ~500,000 single nucleotide polymorphisms (SNPs) using the Affymetrix Mapping 500K SNP chip. 334,923 SNPs passed quality control (QC) procedures29. Validation set I consisted of 808 MS cases and 1,720 controls ascertained through the Partners MS Center in Boston, Massachusetts and genotyped using the Affymetrix 6.0 platform. From this panel, 453 of the 469 SNPs passed QC procedures, and these were used to assess the significant models from the genome-wide screen. Validation set II consists of an independent set of 2,330 MS cases and 2,110 controls ascertained from Brigham and Women’s Hospital in Boston, MA, University of California San Francisco, Washington University, the Accelerated Cure Project out of Massachusetts, Rush University of Chicago, and Cambridge University in the United Kingdom. Eight SNPs were genotyped for this study using the Sequenom MassARRAY iPLEX genotyping platform. Validation set III consists of an independent set of 875 MS cases and 903 controls ascertained as part of the multicenter collaborative GeneMSA study, involving the University of California San Francisco, Vrije Universiteit Medical Center in Amsterdam, University Hospital Basel, and Glaxo SmithKline. This study typed samples using the Illumina Sentrix® HumanHap550 BeadChip genotyping platform6. Because the initial studies used Affymetrix SNP panels, surrogate SNPs were selected from the Illumina platform with close physical proximity and high r-squared value based on the Hapmap CEU data.

Pair-wise linkage disequilibrium (LD) statistics computed for over two million SNPs by the International HapMap Project (posted June, 2006) were used to establish the Caucasian-specific haplotype block boundary for each of the 334,923 SNPs in the IMSGC data set. We defined the boundaries of the haplotype block represented by each IMSGC SNP using an iterative procedure that extends the block boundary sequentially (by SNP) if the D’ measure between the HapMap SNP and the IMSGC SNP is equal to 1. Because the IMSGC SNPs are a subset of all known genomic variants, using HapMap LD statistics in this way provides the larger genomic region (which may harbor susceptibility variants) represented by each IMSGC SNP. 5,137 markers in the IMSGC data set were not represented in the HapMap LD data, and the nearest HapMap marker was used as a surrogate to assess haplotype block boundaries. Marker-gene mappings were generated if a haplotype block that overlaps with any portion of a gene as described by the Ensembl database. IMSGC markers capture 14,236 genes using LD, compared to 13,425 using the markers without accounting for LD.

Using a collection of public data sources that suggest putative gene-gene interaction, we generated a set of roughly 20 million two-SNP models36. These sources include the Kyoto Encyclopedia of Genes and Genomes, the Protein Families database, the Gene Ontology, Reactome, the Database of Interacting Proteins, NetPath, the Genetic Association Database, prior regions of suspected linkage for MS, hand-selected candidate genes, and genes showing differential expression in MS. These resources define gene categories, such as a Gene Ontology term or a KEGG pathway. Because the Gene Ontology is a hierarchical resource with some broad categories, GO terms were used for model generation were restricted to those containing 30 genes or less. For the collection of genes within each category, two-SNP models were exhaustively generated by selecting LD-mapped SNPs (using the LD-Spline procedure above) from two different genes of the category. Models containing two SNPs within the same gene were avoided to prevent the assessment of haplotype effects.

Models were evaluated using conditional logistic regression37 in the trio analysis and logistic regression for case-control analysis. Regression models contained three terms; the additive main effect of each of two SNPs and a multiplicative interaction term. Two test statistics are generated for each model; a model fit statistic (MF) describing the likelihood of the specified model given the data, and a likelihood ratio test statistic (LR) comparing the full model to a reduced model containing only main effect terms. A significant likelihood ratio test indicates that including an interaction term significantly improves the fit of the model. We required both statistics to be significant, further constraining our results set to models with evidence of non-additivity, consistent with Fisher’s description of epistasis38.