Impact of genetic risk loci for multiple sclerosis on expression of proximal genes in patients

Abstract

Despite advancements in genetic studies, it is difficult to understand and characterize the functional relevance of disease-associated genetic variants, especially in the context of a complex multifactorial disease such as multiple sclerosis (MS). As a large proportion of expression quantitative trait loci (eQTLs) are context-specific, we performed RNA-Seq in peripheral blood mononuclear cells from MS patients (n = 145) to identify eQTLs in regions centered on 109 MS risk single nucleotide polymorphisms and 7 associated human leukocyte antigen variants. We identified 77 statistically significant eQTL associations, including pseudogenes and non-coding RNAs. Thirty-eight out of 40 testable eQTL effects were colocalized with the disease association signal. As many eQTLs are tissue specific, we aimed to detail their significance in different cell types. Approximately 70% of the eQTLs were replicated and characterized in at least one major peripheral blood mononuclear cell-derived cell type. Furthermore, 40% of eQTLs were found to be more pronounced in MS patients compared with non-inflammatory neurological diseases patients. In addition, we found two single nucleotide polymorphisms to be significantly associated with the proportions of three different cell types. Mapping to enhancer histone marks and predicted transcription factor binding sites added additional functional evidence for eight eQTL regions. As an example, we found that rs71624119, shared with three other autoimmune diseases and located in a primed enhancer (H3K4me1) with potential binding for STAT transcription factors, significantly associates with ANKRD55 expression. This study provides many novel and validated targets for future functional characterization of MS and other diseases.

Introduction

Multiple sclerosis (MS) is a chronic inflammatory disease of the central nervous system (CNS), which leads to demyelination and neuronal loss. Genetics play a significant role in this disease, the heritability being 0.64 (1). Large efforts have been made to identify genetic features that contribute to the pathogenesis of MS, resulting in a growing list of established disease-associated loci, including human leukocyte antigen (HLA) variants and single nucleotide polymorphisms (SNPs) across the genome (2–4). In addition, some lifestyle exposures affect the susceptibility of MS such as smoking, body mass index, sun-exposure and viral infections (5). Several genetic regions previously known to be associated with MS and other autoimmune or chronic inflammatory diseases, have been fine-mapped with a dense set of SNPs and tested for association with MS within the ImmunoChip Project (3). More than 100 non-HLA SNPs are currently established as MS susceptibility SNPs. Many of these may be causal variants, because fine-mapping compared to conventional genome-wide association studies (GWAS) increases the probability of pinpointing disease-causing variants rather than variants in linkage disequilibrium (LD). Interestingly, there is an overlap of MS-associated variants with other autoimmune diseases ranging from 0.9% for psoriasis to 9.1% for inflammatory bowel disease, Crohn’s disease and primary biliary cirrhosis (3).

Given the non-coding nature of the majority of non-HLA MS SNPs (6), the widespread presence of expression quantitative trait loci (eQTLs) in the genome (7), and the potential impact of gene expression on biological pathways involved in disease, it is likely that a number of MS-associated loci would be involved in regulation of gene expression, acting as eQTLs. A large proportion of the eQTLs present in the human genome seem to be context-specific, i.e. they are active in specific cell types and/or as a response to molecular signaling induced by intrinsic or extrinsic stimuli (6,8,9). Most eQTL studies performed to date, including those focusing on MS variants (10), have used cells from healthy individuals, while studies in the disease context, using patient samples, are still limited. Importantly, in MS chronic inflammation is found in the CNS, with increased numbers of inflammatory cytokines and immune cells in the cerebrospinal fluid, but the inflammatory state is also reflected in the peripheral blood (11) which is more easily accessible.

Therefore, in order to study MS susceptibility variants as eQTLs in their disease and cell-type specific context, we performed whole transcriptome RNA sequencing (RNA-Seq) in peripheral blood mononuclear cells (PBMCs) from 145 patients with MS or clinically isolated syndrome (CIS) that had previously been genotyped on the ImmunoChip SNP microarray (3). Additionally, we included 36 non-inflammatory neurological disease (NINDs) controls, i.e. patients with other neurological diseases in the analysis in order to compare and to identify potential MS-specific effects. Following the eQTL analysis, we characterized the significance of the identified MS eQTLs in lymphoblastic cell lines (LCLs) and primary immune cells (CD4⁺, CD8⁺, B cells and monocytes) sorted from PBMCs obtained from different cohorts of healthy individuals or MS patients. And finally, we added functional evidence for most relevant eQTLs by mapping to enhancer histone marks and predicted transcription factor binding sites. The schematic overview of the approach followed in this study and key results are shown in Figure 1.

Figure 1.

Schematic overview of the study and key results. (A) Immunochip study resulted in 110 non-HLA and 7 HLA markers associated to MS (Supplementary Material, Table S2). (B) We deterimind expression levels in PBMC from 181 patients with RNAseq (Supplementary Material, Table S1) (C) MS-CIS patients eQTL analysis were performed for 109 non-HLA and 5 HLA markers and resulted in 77 significant eQTLs (29 eQTLs associated to HLA markers) (Tables 1 and 2 and Supplementary Material, Tables S3 and S4). (D) Two eQTLs associated to cell propotions were identified: rs2726518 (associated with TET2) was associated with the proportions of monocytes, and rs6881706 (associated with SPEF2) was correlated with the proportions of B cells, CD8+ cells and monocytes (Supplementary Material, Table S6). (E) Colocalization tests performed in 22 MS susceptibility loci. In total, 38 SNP–gene eQTLs were colocalized with the MS locus signal and two SNP–gene eQTLs—rs201202118-XRCC6BP1 and rs533646-AP002954.4.1, were not colocalized (Fig. 2 and Supplementary Material, Figs S2–S5 and Tables S9 and S10). (F) Validation of MS eQTLs in cell-type specific datasets—CD4, CD8 cells (RNA-Seq) and B cells and monocytes (microarray). Of the non-HLA SNP–gene eQTLs found in PBMCs, 40% were significant in at least three additional cell types and 74% were significant in at least one cell type (Figs 3 and 4 and Supplementary Material, Tables S11 and S12). (G) MS eQTL effect changes in monocytes obtained from healthy individuals activated by different stimuli (LPS, 2 and 24 h; IFN-γ, 24 h) compared to unstimulated. eQTL pairs such as rs8042861-IQGAP1 and rs941816-ETV7 showed increased effects in monocytes stimulated with IFN-γ and LPS (2 h), rs2288904-SLC44A2 displayed an increased effect for IFN-γ, and rs2523822 (HLA-A*02)-HLA-H for each of the three stimuli, while rs11052877-CLECL1 showed a decreased eQTL effects after the three stimulations (Fig. 5 and Supplementary Material, Table S13). (H) Disease specificity of eQTL effects—32 out of 73 of the eQTL effects had a slope change <5%, with HLA- B*44 or 45-HCP5 showing the largest (41%) eQTL effect increase. FCRL3 (eQTL gene), FCRL2 and FCRL5 genes in the MS defined locus of rs706015 were significantly downregulated in MS cohort compared to NINDs (Supplementary Material, Fig. S1 and Table S19). (I) Epigenomic mapping of MS-associated variants—Out of the 26 non-HLA eQTLsof MS patients eight overlapped either with an H3K27ac or H3K4me1 active enhancer histone marks or with a DNaseI peak in monocytes, B, CD4⁺ T and CD8⁺ T cells (Fig. 2 and Supplementary Material, Figs S2–S5 and Table S20).

This is, to our knowledge, the first study where the established MS susceptibility variants have been tested for cis-eQTL effects in the disease-specific context by using PBMCs from patients.

Results

Non-HLA and HLA eQTL effects in MS patients

We sequenced the transcriptome of PBMCs from 145 patients with MS or CIS (Cohort 1, Supplementary Material, Table S1), previously genotyped on the ImmunoChip (3), and quantified gene expression using RNA-Seq. Then, we identified pairs of MS SNP–genes by selecting expressed genes within a window of 400 kb on each side of 109 MS-associated non-HLA SNPs (Supplementary Material, Table S2) and 7 MS-associated HLA variants (presence of HLA-DRB1*03: 01, HLA-DRB1*15: 01 and HLA-DRB1*14: 01/14: 02 and absence of HLA-A*02: 01, HLA-B*44 and or 45, HLA-DRB1*07 and HLA-DRB1*01). The window size was selected based on the increased probability of finding cis-regulated genes at closer distances (12,13) and to limit the number of comparisons (Supplementary Material, Tables S3 and S4). For each MS SNP–gene pair identified, we conducted an eQTL analysis. We first regressed out the effects of the relevant confounders [identified in a principal component analysis (PCs), Supplementary Material, Table S5] of gene expression and next, we computed the non-parametric correlation between residuals and the genotype using an additive model. The total number of MS SNP–gene pairs tested was 2223.

We defined eQTL significance using a permutation-based False Discovery Rate (FDR) threshold of 0.01. Table 1 (column A) lists 49 significant non-HLA eQTLs identified in Cohort 1. Table 2 (column A) lists the 28 significant HLA-associated eQTLs from 7 HLA MS risk alleles (Supplementary Material, Table S4 shows the full HLA eQTL analysis results for Cohort 1).

Table 1.

List of MS-associated non-HLA SNPs showing significant eQTL effects

SNP	Chr	Gene	A			B
			Cohort 1			Pooled cohort
			Rho	Bootstrap Rho mean (5%, 95% range)	FDR	Rho	Bootstrap Rho mean (5%, 95% range)	FDR
rs2050568a	1	FCRL3b	−0.29	−0.28 (−0.52, −0.02)	5.6×10−4	−0.26	−0.26 (−0.49, 0.01)	5.2×10−4
rs3748817a	1	FAM213B	0.23	0.23 (−0.04, 0.48)	6.5×10−4	0.22	0.22 (−0.05, 0.47)	3.5×10−3
	1	MMEL1	0.36	0.36 (0.12, 0.57)	8.3×10−6	0.32	0.32 (0.08, 0.53)	1.5×10−5
	1	TNFRSF14	−0.22	−0.22 (−0.46, 0.03)	8.2×10−3	−0.21	−0.21 (−0.45, 0.04)	4.1×10−3
rs17174870	2	MERTK	0.22	0.22 (−0.04, 0.47)	8.6×10−3	0.26	0.25 (0, 0.51)	5.3×10−4
rs7595717a	2	PLEK	0.22	0.22 (−0.04, 0.45)	7.8×10−3	0.27	0.27 (0.02, 0.5)	3.4×10−4
rs842639	2	AHSA2	−0.26	−0.26 (−0.51, 0)	1.7×10−3	−0.30	−0.29 (−0.52, −0.04)	6.4×10−5
rs1920296	3	IQCB1	0.46	0.46 (0.23, 0.67)	<3.6×10−8	0.40	0.39 (0.14, 0.61)	<3.6×10−8
rs2726518	4	TET2c	−0.26	−0.26 (−0.49, −0.01)	1.5×10−3	−0.18	−0.18 (−0.41, 0.08)	0.02
rs7665090a	4	KRT8P46	−0.24	−0.24 (−0.46, 0)	3.3×10−3	−0.28	−0.27 (−0.49, −0.03)	1.9×10−4
	4	MANBAd	−0.24	−0.23 (−0.45, 0.02)e	4.6×10−3	−0.2	−0.2 (−0.43, 0.05)	7.5×10−3
	4	RP11-10L12.6.1d	−0.36	−0.36 (−0.56, −0.13)	8.4×10−6	−0.36	−0.36 (−0.56, −0.12)	1.1×10−6
rs6881706a	5	SPEF2f	−0.22	−0.22 (−0.45, 0.05)	8.8×10−3	−0.22	−0.22 (−0.47, 0.03)	2.6×10−3
rs71624119a	5	ANKRD55	0.30	0.3 (0.05, 0.52)	2.7×10−4	0.30	0.3 (0.06, 0.52)	3.7×10−5
rs11154801	6	AHI1	−0.65	−0.64 (−0.8, −0.43)	<3.6×10−8	−0.64	−0.64 (−0.8, −0.42)	<3.6×10−8
	6	RP3-388E23.2.1d	−0.29	−0.29 (−0.51, −0.05)	4.3×10−4	−0.31	−0.31 (−0.52, −0.07)	2.5×10−5
rs17066096	6	PEX7	−0.25	−0.25 (−0.49, 0)	2.2×10−3	−0.22	−0.22 (−0.47, 0.03)	2.7×10−3
rs941816	6	ETV7	−0.41	−0.41 (−0.61, −0.18)	2.1×10−7	−0.38	−0.38 (−0.59, −0.14)	2.9×10−7
rs706015	7	HOTAIRM1	0.24	0.24 (−0.04, 0.48)e	3.7×10−3	0.20	0.2 (−0.08, 0.45)	7.8×10−3
	7	HOXA3	−0.22	−0.22 (−0.46, 0.02)e	7.6×10−3	−0.17	−0.17 (−0.41, 0.1)	0.03
rs1021156	8	FAM164A	−0.59	−0.58 (−0.74, −0.39)	<3.6×10−8	−0.54	−0.54 (−0.7, −0.32)	<3.6×10−8
	8	IL7d	−0.26	−0.26 (−0.49, −0.01)	1.8×10−3	−0.23	−0.22 (−0.46, 0.04)	2.5×10−3
	8	RP11-578O24.2.1d	−0.51	−0.51 (−0.69, −0.3)	<3.6×10−8	−0.49	−0.49 (−0.67, −0.27)	<3.6×10−8
rs523604	11	PHLDB1	−0.25	−0.25 (−0.49, 0.02)	2.4×10−3	−0.26	−0.26 (−0.49, 0)	5.1×10−4
rs533646	11	AP002954.4.1d,g	0.53	0.53 (0.31, 0.7)	<3.6×10−8	0.52	0.52 (0.31, 0.7)	<3.6×10−8
rs694739a	11	AP003774.1.1d	−0.57	−0.57 (−0.73, −0.37)	<3.6×10−8	−0.58	−0.57 (−0.73, −0.37)	<3.6×10−8
rs11052877	12	CLECL1	−0.26	−0.26 (−0.49, 0)	1.5×10−3	−0.30	−0.3 (−0.53, −0.04)	5.3×10−5
	12	RP11-726G1.1.1d	0.32	0.32 (0.08, 0.53)	9.6×10−5	0.33	0.32 (0.1, 0.53)	8.9×10−6
rs201202118	12	METTL21B	−0.66	−0.66 (−0.79, −0.48)	<3.6×10−8	−0.69	−0.69 (−0.82, −0.52)	<3.6×10−8
	12	XRCC6BP1g	0.3	0.3 (0.05, 0.52)	2.9×10−4	0.31	0.31 (0.07, 0.53)	2.3×10−5
rs4772201	13	CLYBL	0.28	0.28 (0.04, 0.5)e	6.1×10−4	0.22	0.22 (−0.02, 0.45)	2.8×10−3
rs12148050a	14	TRAF3d	0.28	0.28 (0.04, 0.5)	7.9×10−4	0.24	0.24 (−0.01, 0.48)	1.4×10−3
rs8042861	15	IQGAP1	0.24	0.23 (−0.02, 0.48)e	4.2×10−3	0.19	0.18 (−0.08, 0.44)	0.01
rs1886700	16	CDH1	−0.31	−0.31 (−0.52, −0.07)	1.8×10−4	−0.29	−0.29 (−0.51, −0.03)	7.6×10−5
rs7204270	16	MAPK3	−0.35	−0.35 (−0.57, −0.11)	1.6×10−5	−0.34	−0.34 (−0.55, −0.1)	4.3×10−6
	16	RP11-347C12.1.1d	0.33	0.32 (0.07, 0.55)	7.2×10−5	0.33	0.33 (0.08, 0.55)	8.1×10−6
	16	TBX6	0.28	0.28 (0.02, 0.51)	7.1×10−4	0.29	0.29 (0.03, 0.52)	7.3×10−5
	16	YPEL3	0.25	0.25 (−0.02, 0.49)	2.7×10−3	0.25	0.25 (0, 0.49)	7.2×10−4
rs12946510a	17	GSDMB	−0.48	−0.48 (−0.68, −0.24)	<3.6×10−8	−0.48	−0.48 (−0.69, −0.23)	<3.6×10−8
	17	ORMDL3	−0.55	−0.55 (−0.72, −0.34)	<3.6×10−8	−0.54	−0.54 (−0.71, −0.32)	<3.6×10−8
rs4794058	17	AC040934.1d	−0.25	−0.24 (−0.47, 0)	3.1×10−3	−0.3	−0.29 (−0.52, −0.05)	6.0×10−5
	17	NPEPPS	0.23	0.23 (−0.02, 0.47)	5.0×10−3	0.23	0.23 (−0.02, 0.46)	1.9×10−3
	17	TBKBP1	0.53	0.53 (0.3, 0.72)	<3.6×10−8	0.54	0.53 (0.32, 0.71)	<3.6×10−8
rs8070345a	17	RNFT1	−0.24	−0.24 (−0.46, 0.01)	3.3×10−3	−0.24	−0.24 (−0.47, 0.02)	1.4×10−3
rs2288904	19	PDE4A	−0.22	−0.22 (−0.46, 0.03)	7.0×10−3	−0.22	−0.22 (−0.45, 0.04)	3.6×10−3
	19	SLC44A2	0.31	0.31 (0.04, 0.53)	1.9×10−4	0.29	0.28 (0.01, 0.53)	1.0×10−4
rs2283792	22	PPIL2	−0.31	−0.31 (−0.53, −0.07)	1.9×10−4	−0.28	−0.28 (−0.5, −0.03)	1.4×10−4
	22	TOP3B	−0.44	−0.44 (−0.62, −0.24)	3.6×10−8	−0.38	−0.38 (−0.58, −0.15)	2.1×10−7
rs470119	22	CPT1B	−0.31	−0.3 (−0.52, −0.07)	2.1×10−4	−0.28	−0.28 (−0.51, −0.02)	1.6×10−4

Cohort 1, MS specific cohort; Pooled Cohort, Have both MS specific cohort and non-Inflammatory disease (NINDs) cohort.

a

SNP in region implicated to be associated with other autoimmune diseases.

b

Differentially expressed gene, downregulated in MS patients compared to NINDs.

c

After correction for cell-type proportions rs2726518 was found to associate with increased levels of monocytes, but not the TET2 gene.

d

Novel eQTL genes found in PBMC.

e

eQTL genes which are enriched in MS specific cohort.

f

After correction for cell-type proportions rs6881706 was found to associate with increased levels of B cells, CD8⁺ T cells and monocytes, respectively.

g

Not colocalized to disease association signal in the region.

Table 2.

List of MS-associated HLA variants showing significant eQTL effects

HLA allele	Gene	A			B
		Cohort 1			Pooled cohort
		Rho	Bootstrap Rho mean (5%, 95% range)	FDR	Rho	Bootstrap Rho mean (5%, 95% range)	FDR
HLA-A*02:01 (absence)	IFITM4P	−0.50	−0.5 (−0.68, −0.29)	<1.87×10−7	−0.46	−0.45 (−0.65, −0.23)	<1.87×10−7
	HLA-F-AS1	−0.33	−0.33 (−0.55, −0.09)	4.1×10−6	−0.32	−0.32 (−0.54, −0.07)	8.1×10−6
	HLA-H	−0.27	−0.27 (−0.48, −0.03)	1.8×10−4	−0.24	−0.24 (−0.47, 0.01)	9.0×10−4
	MICD	−0.25	−0.25 (−0.47, −0.02)	4.8×10−4	−0.27	−0.26 (−0.49, −0.02)	2.3×10−4
	HLA-K	−0.25	−0.24 (−0.48, 0)	7.0×10−4	−0.25	−0.25 (−0.48, 0)	4.8×10−4
	HCG4a	−0.24	−0.24 (−0.47, 0.01)	1.0×10−3	−0.25	−0.25 (−0.48, 0)	5.2×10−4
	HLA-T	−0.22	−0.22 (−0.47, 0.05)	2.0×10−3	−0.24	−0.23 (−0.47, 0.03)	1.2×10−3
	ZDHHC20P1	−0.25	−0.25 (−0.49, 0)b	5.8×10−4	−0.2	−0.2 (−0.43, 0.05)	6.3×10−3
HLA-DRB1*03:01 (presence)	HLA-DQA1	−0.44	−0.44 (−0.61, −0.24)	<1.87×10−7	−0.47	−0.47 (−0.64, −0.26)	<1.87×10−7
	HLA-DQB1	−0.41	−0.41 (−0.58, −0.2)	<1.87×10−7	−0.43	−0.43 (−0.61, −0.21)	<1.87×10−7
	PPP1R2P1a	0.26	0.25 (−0.01, 0.49)b	4.0×10−4	0.19	0.18 (−0.1, 0.44)	0.01
	HLA-DRB5	−0.25	−0.25 (−0.47, −0.01)	5.1×10−4	−0.26	−0.26 (−0.48, −0.01)	2.6×10−4
	HLA-DQB1-AS1	−0.37	−0.36 (−0.56, −0.14)	3.7×10−7	−0.4	−0.39 (−0.58, −0.17)	1.9×10−7
	HLA-DQB2	0.36	0.35 (0.12, 0.55)	9.4×10−7	0.33	0.33 (0.11, 0.53)	3.7×10−6
	PSMB9	0.22	0.21 (−0.03, 0.44)b	3.1×10−3	0.15	0.15 (−0.11, 0.38)	0.05
HLA-DRB1*15:01 (presence)	HLA-DRB5	0.78	0.78 (0.64, 0.87)	<1.87×10−7	0.77	0.77 (0.62, 0.86)	<1.87×10−7
	HLA-DRB1	0.68	0.68 (0.53, 0.79)	<1.87×10−7	0.66	0.65 (0.49, 0.78)	<1.87×10−7
	HLA-DQB1	0.61	0.61 (0.41, 0.76)	<1.87×10−7	0.56	0.55 (0.34, 0.72)	<1.87×10−7
	HLA-DQB2	−0.51	−0.51 (−0.69, −0.28)	<1.87×10−7	−0.46	−0.45 (−0.65, −0.23)	<1.87×10−7
	HLA-DRB6	−0.48	−0.48 (−0.67, −0.25)	<1.87×10−7	−0.48	−0.47 (−0.67, −0.24)	<1.87×10−7
	HLA-DQB1-AS1	0.46	0.46 (0.24, 0.65)	<1.87×10−7	0.42	0.41 (0.19, 0.61)	<1.87×10−7
	HLA-DQA1	0.39	0.38 (0.14, 0.59)	1.9×10−7	0.34	0.34 (0.1, 0.56)	2.1×10−6
	HLA-DQA2	−0.36	−0.35 (−0.56, −0.11)	9.4×10−7	−0.36	−0.36 (−0.57, −0.13)	5.6×10−7
	TAP2	0.29	0.29 (0.05, 0.52)b	5.0×10−5	0.21	0.2 (−0.04, 0.43)	5.1×10−3
	BTNL2a	0.21	0.21 (−0.05, 0.46)	3.5×10−3	0.19	0.19 (−0.08, 0.44)	8.9×10−3
	NOTCH4	0.22	0.22 (−0.04, 0.45)	2.7×10−3	0.2	0.19 (−0.07, 0.44)	7.8×10−3
HLA- B*44 and or 45 (absence)	HLA-B	0.23	0.23 (−0.02, 0.46)	1.34×10−3	0.26	0.26 (0.01, 0.47)	2.45×10−4
	HCP5a	0.22	0.22 (−0.04, 0.44)b	2.44×10−3	0.13	0.13 (−0.12, 0.38)	0.07

Cohort 1, MS specific cohort; Pooled Cohort, Both MS specific cohort and non-Inflammatory neurological disease (NINDs) cohort.

a

Novel eQTL genes found in PBMC.

b

eQTL genes which are enriched in MS specific cohort.

Filtering eQTL candidate SNPs based on cell proportions

Because the genetic variants themselves may regulate the composition of immune cells within PBMCs (14), the differences in cell-type proportions may explain an observed eQTL effect. To account for that possibility, we first estimated the proportion of different cell types in each sample from the RNA-Seq data using the described method (15). Next, we tested the association of SNPs involved in significant eQTLs (Tables 1 and 2) with the cell proportions (monocytes, NK cells, CD4⁺ and CD8⁺ T cells, and B cells) using Spearman correlation. We identified significant associations to cell proportions (P-value < 0.05) for seven non-HLA SNPs and one HLA variant (Supplementary Material, Table S6). Subsequently, we re-ran the eQTL analysis for these SNPs, while regressing out the cell proportion effects. As a result, two eQTLs were significant associated with cell proportions (Supplementary Material, Table S7); rs2726518 (associated with TET2) was associated with the proportions of monocytes, and rs6881706 (associated with SPEF2) was correlated with the proportions of B cells, CD8⁺ cells and monocytes.

Thus, 47 non-HLA MS SNP–gene pairs and 28 HLA MS SNPs–genes pairs remained significant eQTLs after correcting for cell proportions. Out of these, some eQTL-associated genes, e.g. CPTIB, MANBA, PLEK, METTL21B, AHI1, TNFRSF14, MERTK, IQCB1 and CLECL1, have been reported previously in studies from healthy individuals (8,10,16). Importantly, we also identified eQTLs associated with, to our knowledge, novel genes (mostly non-coding RNAs, antisense transcripts and pseudogenes). Thirty-one out of 47 non-HLA MS-eQTLs (66%) are located within 100 kb of the gene transcription start site from the SNP, although 70.2% of the total eQTLs did not affect the gene closest to the SNP (Supplementary Material, Table S8).

MS eQTLs and linkage disequilibrium

Observed eQTL effects could potentially be caused by SNPs in linkage disequilibrium (LD) with the established MS-associated SNPs. Using fine-mapped genotypic information from the ImmunoChip study, we investigated all SNPs with r²>0.5 (referred to as LD SNPs) in relation to the established top-associated MS SNPs showing eQTL effects in our data. Many of the LD SNPs showed similar or more significant eQTL effects on the same genes as the established MS SNPs (Supplementary Material, Table S9). About 49 out of 71 SNP–gene (including HLA) loci were found to have at least one SNP with better eQTL P-value compared to that of MS SNP–gene. To identify an appropriate candidate eQTL SNP in each locus, we first performed a conditional analysis where, for each MS SNP and for each associated LD SNP, the effect of each associated LD SNP was first regressed out and then residuals were correlated with the MS SNP. None of the MS SNP eQTL effects remained significant when conditioned on any of the SNPs with r²>0.5. Subsequently, we conducted the reverse analysis for each eQTL. For each MS SNP and for each associated LD SNP, we first regressed out the eQTL effect of the MS SNP and then residuals were correlated with the LD SNP as described previously.

Five regions contained LD SNPs with significantly associated (P < 0.05) eQTL effects (Fig. 2 and Supplementary Material, Figs S2–S5). These LD SNPs were associated with the expression of ORMDL3/GSDMB, TNFRSF14, METTL21B/XRCC6BP1, FCRL3 and HOTAIRM1 suggesting additional independent eQTL effects. With the exception of FCRL3, all these genes were located in regions that had been fine-mapped on the ImmunoChip.

Figure 2.

MS associations, LD and mapping to DNaseI peaks, active enhancer chromatin marks and TFBS of ImmunoChip SNPs in the METTL21B/XRCC6BP1 region. The ImmunoChip association analysis results (red dots) are shown for each of the SNPs. SNPs labelled in red have transcription factor binding site based on RegulomeDB—STAT1 (rs10783848), NTAF (rs1875124), BCL6B (rs11172344), SRF (rs11172343), IRF3 (rs923829), HELIOSA (rs10877015), FOXP (rs12423195), SOX4 (rs12423195), FOXL1 (rs10877012) and E2F2 (rs4646536). The top panel shows regions with H3K4me1 and H3K27ac chromatin mark peaks (red) and DNaseI peaks (black) in the different cell types analyzed in this study (PBMCs, Monocytes, CD8+ and CD4+ T cells and B cells). SNP–gene correlation P-value for two genes—METTL21B and XRCC6BP are marked as solid lines. XRCC6BP eQTL association signal is not colocalized with the SNP disease association signals.

For the remaining 21 LD regions included in the conditional analysis, none of the LD SNP eQTL effects remained significant when conditioned on the MS SNP (Supplementary Material, Table S9). We conclude that with few exceptions our results are consistent with ImmunoChip top-associated MS SNPs affecting MS risk by controlling expression levels. However, with an expression cohort with moderate sample size, we cannot exclude the possibility that the eQTL effect we observed for the MS-associated SNPs is due to SNPs in LD.

MS eQTLs and colocalization with MS risk

For the refinement of causal variants based on disease-associated genetic variants and gene expression, we performed colocalization tests in 22 MS susceptibility loci including HLA loci (the number was obtained after applying the criteria described in Materials and Methods), corresponding to 40 eQTL-gene pairs. The colocalization method ‘COLOC’ provides an approximate Bayes factor to estimate a posterior probabilities that a SNP is causal in both immunochip disease association and eQTL analysis. Thirty-eight SNP–gene eQTLs were colocalized with the MS locus signal and two SNP–gene eQTLs—rs201202118-XRCC6BP1 (Fig. 2) and rs533646-AP002954.4.1, were not colocalized with disease association signal from ImmunoChip study (Supplementary Material, Table S10).

Validation and exploration of MS eQTLs in cell-type specific datasets

Considering (a) the context-specific nature of eQTLs, and (b) that most eQTLs are expected to be significant in healthy individuals as well, we aimed to characterize significant eQTLs activity in different cell types using published and novel datasets. We characterized the identified eQTLs in lymphoblastic cell lines (LCLs), naïve B cells, naïve monocytes, as well as monocytes stimulated with LPS (2 and 24 h) or IFN-γ (24 h) sampled from healthy volunteers for which we obtained genotypes and expression data from public datasets (7,8,16) (Supplementary Material, Table S11). We also used a new small cohort consisting of sorted CD4⁺ and CD8⁺ T cells from MS patients and healthy controls (Table 3).

Table 3.

Summary of different datasets used for genotying and screening phase (PBMC) and validation phase (PBMC-derived immune cell types)

Cohort	Tissue	MS	NINDs	HC	Expression type (source)	Analysis results and data
Pooled Cohort [Cohort 1 (MS)+Cohort 2 (NINDs)]	PBMC	145	36a	0	RNA-Seq	Tables 1 and 2 and Supplementary Material, Tables S1, S3, S4 and S15–S18
Cell-type specific	CD4a	19	0	8	RNA-Seq	Figures 3 and 4 and Supplementary Material, Tables S12, S13 and S23
Cell-type specific	CD8a	17	0	7	RNA-Seq	Figures 3 and 4 and Supplementary Material, Tables S12, S13 and S23
Cell-type specific	B-cells	0	0	233	Microarray (16)	Figures 3 and 4 and Supplementary Material, Tables S12 and 13
Cell-type specific	Monocytes	0	0	414	Microarray (8)	Figures 3–5 and Supplementary Material, Tables S12–16
Cell-type specific (stimulation)	Monocytes (IFN gamma)	0	0	367	Microarray (8)	Figures 3–5 and Supplementary Material, Tables S11–S13
Cell-type specific (stimulation-time point)	Monocytes (LPS—2 h)	0	0	261	Microarray (8)	Figures 3–5 and Supplementary Material, Tables S11–S13
Cell-type specific (stimulation-time point)	Monocytes (LPS—24 h)	0	0	322	Microarray (8)	Figures 3–5 and Supplementary Material, Tables S11–S13
Cell-type specific (1000 genome project)	EBV-transformed LCLs (B cells)	0	0	234	RNA-Seq (7)	Figures 3 and 4 and Supplementary Material, Tables S11 and S12
Additional PBMC dataset (MS and NINDs)	PBMC	23	7	0	Human Genome U133 plus 2.0 arrays (18)	Supplementary Material, Table S19
Global Cohort Association study (Cases and Controls)	Blood (DNA)	17465		30385	ImmunoChip custom genotyping array (3)	Supplementary Material, Table S10
Joint Analysis cohort (Cases and Controls)	Blood (DNA)	29300		50794	ImmunoChip custom genotyping array (3)	Supplementary Material, Table S10

LCLs, lymphoblastoid cell lines.

a

Sample size estimation and power analysis done for these cohorts.

We performed eQTL analyses with the same methodology described earlier. Gene expression in the LCLs and primary T cells had been analyzed using RNA-Seq, which allowed us to detect and validate several of the non-coding RNA transcripts or transcribed pseudogenes found in our initial analysis in the PBMCs. The gene expression data from sorted primary immune cells from healthy individuals had been obtained through microarray analyses, where probes were either lacking or expressed below detection limit for some of the eQTL-affected genes. For confirmation of an eQTL effect on a given gene-locus pair, we accepted a nominal P-value of <0.05 if the effect shows the same direction as observed in Cohort 1. More than 40% of the non-HLA SNP–gene eQTLs found in PBMCs was identified significant in at least three additional cell types. Furthermore, 74% were identified significant in at least one cell type (P-value 3.55×10⁻⁷). Some eQTLs were general for all cell types, while most of them were specific for one or a few cell types or conditions (Figs 3 and 4 and Supplementary Material, Tables S12–S16). Importantly, we acknowledge that the analysis using the T cell cohort had low statistical power due to its size; however, the number of replicated eQTLs (P-value < 0.05) is higher or equal to the expected numbers from a power analysis.

Figure 3.

Graphical overview of MS-associated non-HLA eQTLs in immune cell types. Validation and exploration of MS non-HLA eQTLs in LCLs, primary B cells, primary monocytes [unstimulated, stimulated with IFN-γ (24 h), LPS (2 h), LPS (24 h)], from healthy individuals, as well as CD4⁺ and CD8⁺ T cells from MS patients and healthy controls. A nominal P-value of <0.05 was accepted as validation of a significant eQTL in PBMCs. *A proxy SNP was used in the analysis of primary B cells and monocytes (listed in Supplementary Material, Table S24). eQTL SNPs (rs2050568, rs2726518 rs12946510 and rs201202118) were excluded as original snps or tag snps were not obtained for the corresponding studies. +/− indicates the gene expression increase/decrease with respect to the risk allele.

Figure 4.

Graphical overview of MS-associated HLA eQTLs in immune cell types. Validation and exploration of MS HLA eQTLs in LCLs, primary B cells, primary monocytes [unstimulated, stimulated with IFN-γ (24 h), LPS (2 h), LPS (24 h)], from healthy individuals, as well as CD4+ and CD8+ T cells from MS patients and healthy controls. *rsids used as tag SNPs for the HLA alleles in LCLs, B cells and monocytes. +/− indicates the gene expression increase/decrease with respect to the risk allele.

In the analysis of eQTL effects in the HLA region, we found that the HLA-DRB1*15:01 allele or its tag SNP, rs3135388, significantly influenced expression of the DRB1, DRB5 and DQB1 genes throughout the tested datasets except T cells (Fig. 4 and Supplementary Material, Table S12).

Influence of LPS and IFN-γ stimulation on MS eQTLs in monocytes

To further characterize the context-specific nature of the eQTLs in monocytes, which are thought to be involved in MS pathogenesis (17), we analyzed the behavior of the MS eQTLs in monocytes obtained from healthy individuals (8) activated by different stimuli (LPS, 2 and 24 h; IFN-γ, 24 h) compared to unstimulated, by comparing the changes in effects, measured as rho-values (Fig. 5 and Supplementary Material, Table S13). Most increases in MS eQTL effects (rho) were observed upon IFN-γ stimulation, while a certain predominance of decreased eQTL effects was seen upon stimulations with LPS as compared with unstimulated cells. In accordance with the pattern observed in PBMCs from MS patients, eQTL pairs such as rs8042861-IQGAP1 and rs941816-ETV7 showed increased effects in monocytes stimulated with IFN-γ and LPS (2 h), rs2288904-SLC44A2 displayed an increased effect for IFN-γ, and rs2523822 (HLA-A*02)-HLA-H for each of the three stimuli, while rs11052877-CLECL1 showed a decreased eQTL effects after the three stimulations.

Figure 5.

Heat map of changes in MS eQTL effects upon stimulation of monocytes. The changes are measured as rho values and are depicted in percentage relative to the unstimulated (naïve) condition. Only SNP–gene pairs with nominally significant (P<0.05) effects in the naïve state are included in the heatmap. Graphs are llustrating eQTL effect size (beta-value) changes for the different stimuli in monocytes, showing examples of expression of four genes: ETV7, SLC44A2, IQGAP1 and FAM164A.

Disease-specific genes in MS, CIS and NINDs

To understand the role of MS-SNP regulated genes in the disease, gene expression levels of MS were compared to those of NINDs and CIS (Supplementary Material, Table S14). Thirty-nine out of 45 differentially expressed genes were downregulated in MS compared to NINDs (adjusted P-value < 0.05) and no genes were significantly differentially expressed between MS and CIS samples (Supplementary Material, Table S14). FCRL3 (eQTL gene), FCRL2 and FCRL5 genes in the MS defined locus of rs706015 were significantly downregulated in MS cohort compared to NINDs (Supplementary Material, Fig. S5).

Disease specificity of eQTL effects

While we only observe one eQTL gene, FCRL3, to be differentially expressed, we aimed to investigate if the eQTL effects are disease specific by displaying an increased association in MS, or rather general and independent of the disease status.

To this end, we pooled RNA-Seq and genotyping data from 36 NINDs (Cohort 2) and MS (Cohort 1) and regressed out the effect of disease type (MS/NINDs) in the gene expression profiles. Then, we compared results from this Pooled cohort with those from Cohort 1 [Tables 1 and 2 (panel B); the complete analysis results for the Pooled Cohort are found in Supplementary Material, Tables S15 and S16]. All non-HLA and HLA eQTLs except two (HLA-DRB1*03:01-PSMB9 and HLA-B*44/45–HCP5) remained significant in the Pooled Cohort. We calculated the percentage of relative increase or decrease of eQTL effects in the MS group (Cohort 1) as compared to when patients with all diagnoses are included (Pooled Cohort) (Supplementary Material, Fig. S1A and B and Tables S17 and S18). The results indicated that 32 out of 73 of the eQTL effects had a slope change <5%, with HLA- B*44 or 45-HCP5 (Supplementary Material, Fig. S1B) showing the largest (41%) eQTL effect increase. We lack power to quantify accurately the change of the slope per eQTL, however, we observed that for a set of HLA eQTLs the association strength (absolute value of rho) is reduced when adding non-MS individuals (Supplementary Material, Table S17) despite the increase of the number of individuals (P-value < 0.05). In non-HLA eQTLs, we also observed for many eQTLs a reduction of association strength, but it was not statistically significant. However, to further test the validity of increase or decrease of eQTL effects in MS compared to other patients, we computed eQTLs in a small previously analyzed cohort with MS and NINDs (OND) patient samples taken from the same study population as the Pooled Cohort (18) and we found that the change of the effect between MS and NINDs was in the same direction for a significant set of eQTLs (P-value <0.05) (Supplementary Material, Table S19).

In order to address if our results were affected by realapse–remission status, we investigated the eQTLs reported in Tables 1 and 2, stratifing the MS and CIS patients based on realapse–remission status. Among patients in remission (n = 101) three of the eQTLs (HLA-DRB1*03:01-PSMB9, rs11154801- RP3–388E23.2.1 and rs8070345- RNFT1) were not significant while all the others remained significant. However, these three eQTLs were all significant among patients in remission (n = 29), indicating that they may be specific to situation of active inflammatory response (data not shown).

Chromatin marks

To functionally characterize the MS eQTLs identified, we performed a bioinformatics screening of non-HLA MS eQTL loci using chromatin marks and DNaseI peaks, as well as analysis of predicted TFBS mapping from ENCODE data (19) and Roadmap Epigenomics data (20) from healthy individuals. Out of the 26 non-HLA eQTLs in PBMCs of MS patients eight overlapped either with an H3K27ac or H3K4me1 active enhancer histone marks or with a DNaseI peak in monocytes, B, CD4⁺ T and CD8⁺ T cells. Two non-HLA eQTLs also overlapped with a predicted TFBS for PIT-1 and STAT1, while NKX2–3, NKX2–8, Irx5, Irx3 and GRHL1 were predicted to bind SNPs without any other overlapping marks (Supplementary Material, Table S20).

Discussion

We have performed an eQTL analysis of the top-associated MS risk loci (3) by studying cis-eQTL effects in 800-kb windows in PBMCs from patients with MS. We found 35 MS loci that displayed eQTL effects influencing one or more genes within the analyzed regions. A majority of the eQTL effects remained significant when adding data from non-MS patients, however, >70% of the effect sizes changed when analyzed in only MS/CIS patients, with a predominance of increased effects in MS cohort for HLA-related eQTLs (Supplementary Material, Fig. S1).

We replicated the eQTLs identified from PBMCs in datasets from LCLs, primary B cells and monocytes from healthy individuals, including monocytes activated by LPS and IFN-γ as well as in CD4⁺ and CD8⁺ T cells from MS patients and healthy controls. In total, 74% of the non-HLA eQTLs were replicated in at least one cell type. We replicated known observations from healthy individuals (8,10,16) as well as reported previously unidentified eQTLs, especially for HLA genes, long non-coding or anti-sense RNA genes and pseudogenes. Interpretation of comparisons with previously published eQTL studies should however be made with caution, considering the differences in the SNPs or proxy SNPs, expression platforms and number of study subjects included. Importantly, considering the overlap in risk genes between autoimmune diseases, we consider that the results presented here are relevant for other autoimmune diseases (Table 1).

Roles of MS eQTLs in different immune cell types

B cells and monocytes are found in MS lesions and their roles in the disease process are being investigated (21–24), however, there are plenty of studies pointing toward an important role of T cells in the disease (2,3). Therefore, we presented RNA-Seq based eQTL data from a small cohort of CD4⁺ and CD8⁺ T cells, which allowed us to validate a few of the eQTL findings from PBMCs. Despite the small sample size, two of the eQTL effects found in CD4⁺ and CD8⁺ T cells (rs7204270-YPEL3, rs8070345-RNFT1) could not be replicated in any other cell type, thus suggesting T cell specificity (Fig. 3 and Supplementary Material, Table S11). Larger sample sizes will be needed for a more robust exploration of MS SNPs as eQTLs in different T cell subsets.

After integrating chromatin marks and TFBS data, and taking the different validation strategies together, we found that the AHSA2 eQTL, rs842639, had perfectly matching results for DNase I and histone marks throughout all tested datasets, confirming its eQTL effect. For the remaining MS eQTLs histone marks, DNaseI and TFBS data overlapped with eQTL results in the different tissues to a varying extent.

To describe some examples, rs3748817 was found as an eQTL for FAM213B and MMEL1 in PBMCs from MS patients. However, the behavior of this eQTL effect under different stimulatory conditions of monocytes varied (no eQTL for FAM213B in 24-h stimulation with LPS), and it also changed in different PBMC-derived cell types (MMEL1 has no eQTL effect in B cells, and no eQTL effect for any of the genes in T cells). An active enhancer histone mark was also found for this SNP in monocytes (Supplementary Material, Table S20 and Fig. 4). rs523604, an eQTL for PHLDB1 in MS PBMCs, was found in healthy monocytes, however not supported by any other data for this cell type, while it mapped to active enhancer in healthy PBMCs, B cells and CD8⁺ T cells. We confirmed the eQTL for CLYBL, rs4772201, in healthy B cells, while active enhancer histone marks and DNase I peaks were found throughout the analyzed cell types (Supplementary Material, Fig. S2 and Table S20). Most of the TFs with predicted binding to eQTL SNPs have described functions during development (mainly neuronal) and cell cycle. Only STAT1 (predicted to bind to rs71624119, eQTL for ANKRD55) and NKX2–3 (predicted to bind rs11154801, eQTL for AHI1 and RP3–388E23.2.1) have known functions in immunity.

Detection of shared causal variants

A question of importance is whether the disease-associated SNP in a region is the one driving an observed eQTL effect and therefore possibly responsible for disease risk. We show that for at least 6 of the 24 testable MS SNPs there are other SNPs in the region with stronger eQTL effect in PBMCs from MS patients (Supplementary Material, Table S9). The most significant disease-associated SNPs in the global analyses in these regions may simply have other functional effects than eQTLs in PBMCs, or alternatively, previous association analysis and fine mapping has failed to pinpoint the correct SNP for MS susceptibility (Fig. 2 and Supplementary Material, Figs S2–S5). Additionally, for every gene–SNP pair, we conducted a colocalization analysis in a window span of 800 kb around risk SNP, to test whether the eQTL variants colocalizes with the disease association variant. Our study suggest that 95% of the testable disease loci were colocalized to the eQTL gene signal (Supplementary Material, Table S10). Larger cohorts that enable a more thorough dissection of the LD structure and regions with even coverage of SNPs will be needed in order to conclude exactly which SNPs are necessary for the effect, although it is likely that haplotypes, rather than single SNPs are the drivers (25).

Modification of eQTL effects by context

Current data from eQTL analyses in different immune cell types show that a large proportion of eQTL effects are highly context-specific, i.e. they vary depending on the signaling pathways activated by different kinds of stimuli, such as cytokines affecting the particular cell type (6,9,10). MS is an inflammatory disease and patients present with perturbed peripheral blood cytokine levels (26). It is thus likely that the differences in eQTL effect sizes in MS patients as compared to the full cohort with NINDs included (Supplementary Material, Fig. S1) are at least in part due to the exposure of PBMCs to inflammatory cytokines in MS patients. One of the gene with the highest eQTL effect increase in MS/CIS patients, was PSMB9, a proteasome subunit involved in processing of class I MHC peptides (27) and TAP2 which help in antigen transport and function as a molecular scaffold for the final stage of MHC class I folding (28).

We analyzed the MS eQTLs in monocytes from healthy individuals stimulated with LPS and IFN-γ, which may to a degree mimic an inflammatory context similar to MS (29). Both are potent stimuli that boost cytokine production at 24 hours, thus activating additional signaling pathways in the monocytes. We have used the stimulation-specific eQTL data as a potential approximation of the context of chronic inflammation in MS patients, and for some of the genes with increased eQTL effects in MS, such as IQGAP1 or CLECL1, we could indeed observe a similar pattern upon stimulation of monocytes with LPS and/or IFN-γ.

Gene functions and potential relevance in MS disease process

The non-HLA MS genes eQTLs in PBMCs represent a wide variety of functions and processes. Several are involved in the immune system, such as the well-established IL7, TRAF3 (30,31), MERTK (32), CLECL1 (33) and TBKBP1 (34,35). A large proportion have unknown or poorly studied functions (e.g. FAM164A, FAM213B) and another group of genes are the non-translated ones, such as anti-sense RNA genes (RP3–388E23.2.1, AP002954.4.1, CTD-2320O4.2.1, AP003774.1.1) and pseudogenes (RP11–578O24.2.1, RP11–726G1.1.1, AC040934.1, KRT8P46, RP11–10L12.6.1). Even though pseudogenes do not produce functional proteins, there is increasing evidence for regulatory roles for transcribed pseudogenes (36–38). In this study, we have observed significant downregulation of FCRL3 in MS compared to NINDs. SNPs in the FCRL3 locus have been identified as risk SNPs for MS and other autoimmune diseases such as autoimmune thyroid disease (Supplementary Material, Table S2). We also identify FCRL3 as one of the eQTL genes in this study. Furthermore, there is an experimentally validated interaction between FCRL3 and PTPN22 (risk gene for rheumatoid arthritis) identified by the Protein–protein interaction (PPI) network provided in STRING database (39). Therefore, as it is differentially expressed it is a relevant candidate to be further characterized at the protein level. Recently, the possible role of ANKRD55 in neuroinflammation was studied at protein expression level in neuron and microglia cultures based on an animal model of experimental autoimmune encephalomyelitis (EAE) (40). In our study, the eQTL SNP for ANKRD55 is found to be located in an active regulatory region (H3K4me1) in PBMC and four PBMC-derived immune cell types with a predicted score for binding of STAT1 transcription factor and is also found to be associated to other autoimmune diseases such as Crohn's disease, Rheumatoid Arthritis ad Juvenile Idiopathic Arthritis (3,41). Finally, AHI1 (42,43), MANBA and MMEL1 (44,45) have potential roles in the CNS, but could also be important in other cell types. Some genetic variants regulate the expression of several genes, however, the effects observed for the different genes may vary largely depending on the cell type that is being studied, and therefore it is not clear which gene drives the disease association.

eQTL effects in the HLA-DRB1*15:01-containing haplotypes have been investigated previously (46–49). We replicated the differences in expression of HLA-DRB1, DRB5, DQB1 (46,48) and DQA1 (46) in relation to HLA-DRB1*15:01, but we also find novel associations (e.g. HCP5 and HLA-B association to HLA-B*44 or 45). As HLA-DRB5 is present exclusively on HLA-DRB1*15 and HLA-DRB1*16-containing haplotypes, its expression has served as a positive control for the HLA-DRB1*15 allele throughout all the tested cell types (Fig. 4). Presence of HLA-DRB1*15:01 was significantly correlated with increased expression of DRB1 in PBMCs, B cells and monocytes. Importantly, we also found genes affected by the presence of the DRB1*03:01 risk allele (HLA-DQA1, HLA-DQB1, PPP1R2P1, HLA-DRB5, HLA-DQB1-AS1, HLA-DQB2, PSMB9) and the A*02:01 protective allele (IFITM4P, HLA-F-AS1, HLA-H, MICD, HLA-K, HCG4, HLA-T, ZDHHC20P1), which to our knowledge, have not been reported before (Fig. 4).

Extensive efforts are presently being made to map eQTLs in different healthy immune cell subsets and stimulatory contexts. Disease-associated eQTLs in a disease-specific context is an important additional aspect that should be considered in light of our results on increased strength effect (absolute rho) in HLA-related eQTLs. By using MS patient samples, we provide an initial step of insight into potential SNP–gene relationships that lay a basis for a more extensive mapping in the context of MS. We present a detailed overview through which we provide a list of putative MS genes, confirm that these genes are not necessarily those located closest to the SNPs, vary in their response depending on the cell subset examined, and are influenced by disease or stimulatory conditions. This study serves as a resource of information for further functional exploration of genes and signaling pathways that are affected by MS-associated polymorphisms.

Materials and Methods

We have tested correlation between genotype and expression levels in several different cohorts which are briefly described in Table 3.

Patient selection in screening phase for PBMCs

For the screening-phase of the project, we obtained RNA and DNA from PBMC of 181 subjects collected between 2001 and 2010, at the Neurology Clinic of the Karolinska University Hospital, Solna, Sweden. The cohort consisted of 117 MS patients fulfilling the 2001 McDonald criteria (for consistency later revisions to these guidelines were not considered) (50), 28 patients with clinically isolated syndrome (CIS) and 36 patients were NINDs, i.e. patients suffering from neurological diseases without an inflammatory state, such as neuralgia, paresthesia, sensory symptoms, vertigo, tension headache (Supplementary Material, Table S1). About 24 out of 28 CIS patients were later diagnosed to have MS.

Preparation of PBMCs, CD4⁺ and CD8⁺ T cells

We isolated PBMCs immediately from blood samples taken with sodium citrate-containing cell preparation tubes (BD Vacutainer™ CPT™ Tube, Becton-Dickinson, Franklin Lakes, NJ), in the screening phase, or using a standard Ficoll (GE Healthcare, Little Chalfont, UK) procedure, in the validation phase. We separated the cells with density gradient centrifugation, collected them from the interphase, washed twice in Dubelcco’s phosphate buffered saline (PBS) and either froze as a pellet at −80° C for subsequent RNA extraction, or prepared the cells for sorting. Sorting of the CD4⁺ and CD8⁺ T cell population was performed by adding fluorochrome conjugated antibodies against human CD4 (clone SK3, APC-conjugated, Becton Dickinson), CD8 (clone SK1, FITC-conjugated, Becton Dickinson) and CD3 (clone UCHT1, PE-conjugated, BD Bioscience) to the negative fraction obtained after manual sorting for CD14⁺ monocytes and using a MoFlo high speed cell sorter.

RNA isolation from PBMCs, CD4⁺ and CD8⁺ T cells

We performed extractions of PBMC RNA to be used for sequencing using the RNeasy Mini Kit (Qiagen, Hilden, Germany), and included DNase digestion, following the manufacturer’s instructions. We eluted the RNA in 10 mM Tris–HCl, 1 mM EDTA buffer, pH 7.0 (Ambion, Austin, TX) and stored it in −80° C. We measured the RNA concentration on a NanoDrop ND-1000 Spectrophotometer (NanoDrop Technologies, Wilmington, DE), RNA was extracted from CD4⁺ and CD8⁺ T cells with Trizol (Invitrogen) using the miRNeasy Kit (Quiagen) according to manufactureŕs recommendations, and analyzed the RNA integrity with Bioanalyzer (Agilent Technologies, Santa Clara, CA). All samples had an RNA integrity number (RIN) >8.

DNA extraction and genotyping

For DNA extraction, we collected blood in EDTA tubes (Becton-Dickinson), which were stored at −20° C for subsequent procedures. We extracted DNA from whole blood samples using QIAamp DNA Blood Maxi Kit (Qiagen, Hilden, Germany), according to the manufacturer’s protocol, and stored it at −20° C. An Illumina Infinium custom array, designed specifically to fine-map regions with established autoimmune-associated risk loci, as well as to test overlap in association between different autoimmune diseases, was used for genotyping. Genotyping of DNA samples in the screening phase (as a part of the ImmunoChip Project) was performed at the Wellcome Trust Sanger Institute and at the University of Miami (UM), Miller School of Medicine (described in detail in 3). Of relevance for the study described in this paper, we included the 109 established non-MHC risk loci typed on the ImmunoChip and showing the strongest independent association to MS in each region (3) (Supplementary Material, Table S2), in the subsequent eQTL analysis.

Samples used for validation and exploration of MS eQTLs in CD4⁺ and CD8⁺ T cells were genotyped with the Illumina custom array called MS replication chip at Miller School of Medicine, UM as part of an IMSGC project (51).

Imputed HLA types

For the samples from PBMC, we imputed HLA types with HLA*IMP:02 (52), using genotypes obtained as described earlier (ImmunoChip for cohort 1 and 2 and MS replication chip for CD4⁺ and CD8⁺ T cells). Alignment and phasing was done against the HapMap CEU panel. All samples included in subsequent analyses had an imputation quality score >0.7 for the HLA alleles of interest. Both the ImmunoChip and MS replication chip has extensive coverage of SNPs in the MHC region.

Transcriptome analysis in PBMC

cDNA libraries for PBMC samples (Cohort 1 and Cohort 2) were prepared using Illumina TruSeq kit (Illumina, San Diego) and sequenced on an Illumina HiSeq 2000 machine. A total of 100 base long paired-end sequence reads were generated on an average sequence depth of 36 million reads per sample for 181 samples (Table 3). The reads were mapped to the H. sapiens reference genome (NCBI v37, hg19) using STAR aligner (53) and then HTSeq tool (54) was used to quantify counts per gene, applying the default parameters in each case. Conditional Quantile Normalization (CQN) method (55) was used to normalize the gene count datasets accounting for the GC content and length biases.

Transcriptome analysis in CD4⁺ and CD8⁺ T cells

Of total RNA, 500 ng with RIN >8.0 was subjected to the Illumina TruSeq Stranded mRNA Library Preparation Protocol with Dual Indexes (Cat. No: RS-122–2013). Libraries were quantified using the Kapa Library Quantification Kit (Cat. No: KK4824). Sequencing was carried out on the Illumina Hiseq 2500 to generate 75 bp paired-end data with an average of 10M reads above Q30. The sequence reads were mapped to hg19 reference with Tophat and HTseq was use to quantify counts-per-gene. The expression data from CD4⁺ cells (23 samples) and CD8⁺ cells (26 samples) was normalized with CQN method as described in PBMC analysis. Sample size estimation and power analysis was done for both the expression dataset using ‘pwr.r.test’ function implemented in ‘pwr’ R package.

Expression and genotype data from LCLs

We downloaded genotypes from the 1000 Genomes Phase 1 dataset from http://www.1000genomes.org. A detailed description of the dataset is provided in the study (56) (Table 3). We used tagging SNPs for the three HLA variants to be validated in our analysis (all were estimated in Swedish population; 979 individuals with HLA*IMP: 01 imputed genotypes): rs2970533 A tagging HLA-A*02:01 (r²= 0.83), rs2596551 C tagging HLA-B*44 and/or 45 (r²=0.89), rs9271366 G tagging HLA-DRB1*15:01 (r²=0.98) and rs2187668 T tagging HLA-DRB1*03:01 (r²=0.96). rs201202118 was not present in the 1000 Genomes dataset and rs10431552 (r²=1) was used as a proxy.

We downloaded RNA-Seq files in FASTQ format for EBV transformed LCLs from individuals included in the 1000 Genomes Phase 1 dataset http://www.ebi.ac.uk/arrayexpress/experiments/E-GEUV-1/samples/. More details about this public RNA-Seq dataset are found in (7). The reads were mapped to the H. sapiens reference genome (NCBI v37, hg19) using STAR aligner (53) and normalized using the CQN method. In total, 234 individuals (105 females and 129 males) with genotype and expression data were included in this study.

Expression and genotype data from B cells and monocytes

Both microarray-based mRNA expression and genotype data were obtained from the Fairfax et al. studies (8,16) for B cells (233 samples), unstimulated monocytes (414 samples) and stimulated monocytes—IFN-ɣ (367 samples), LPS 2 h (261 samples) and LPS 24 h (322 samples). All the samples collected in this study were obtained from healthy individuals (Table 3). The microarray expression data was obtained as log-2 transformed normalized data. For the missing MS-associated SNPs in the Fairfax datasets, proxy SNPs with high LD (r²>0.8) were identified using the SNAP proxy search tool (57) and were used for the eQTL analysis (Supplementary Material, Tables S12 and S13). rs22837922, rs7665090 and rs4701192 were not present on the genotyping chip used in the study, and no appropriate proxy SNPs were available.

Expression and genotype data from PBMC array

Expression data for PBMC was obtained from the study (18) and 30 samples were genotyped as part of an IMSGC project (51). The expression data was reanalyzed using the latest array annotation data for HGU133Plus2 and expressed probes (after background correction) were reported in RMA measures after quantile normalization.

eQTL analysis of MS-associated SNPs

Different approaches implemented toward the transcriptome analysis for various cohorts used in this study (Table 3) are described in Supplementary Material. For each gene–SNP association, we did a two-level step analysis to compute the strength and significance of the association. First, we regressed out the effects of a selected set of variables based on the result of component-based analysis and, depending on the dataset, available meta-information (Supplementary Material, Table S5). Secondly, residuals were correlated (non-parametric Spearman correlation) with genotype information considering a genetic additive model. For each SNP analyzed, we studied gene–SNP associations for the genes that were located 400-kb upstream and downstream of the SNP. The size of the genomic window was chosen based on the increased probability of finding cis-regulated genes at closer distances up to 400 kb (12,13) and with the aim to limit the number of comparisons. Only genes that were expressed in at least 85% of all the samples were included in the analysis.

To estimate the significance of each gene–SNP association, we computed a permutation-based P-value, where 15 000 permutations were done over the residuals of the regress-out model in the CQN transformed expression dataset. False discovery rate (FDR) was computed using a non-parametric method described elsewhere (58); in all cases monotonicity was enforced and FDR was capped at 1.

Gene–SNP associations with FDR < 0.01 were considered significant. Importantly, we selected a conservative selection criteria aiming for high specificity and low sensitivity in order to minimize false positives. In this study, the eQTL analysis was done separately for HLA region and non-HLA region. For a given SNP or allele, the minimum number of samples carrying homozygous minor allele to implement the eQTL analysis was limited to three. HLA alleles which were reported as risk for MS in immunochip study, such as absence of HLA-B*44:02, absence of HLA-B*38: 01 and absence of HLA-B*55:01 were removed from the analysis due to low number of samples in one particular allele group.

For the regression model in the MS datasets (PBMC, CD4⁺ and CD8⁺ samples), we selected the variables based on Principal Component Analysis (PCA) from the meta-information available for each dataset (batch of RNA-Seq library preparation, age of patient, gender, disease-type, clinical course of MS and CIS and Interferon treatment). In the MS PBMC dataset (Cohort 1) and NINDs dataset (Cohort 2), we selected batch of RNA-Seq library preparation as covariate (Supplementary Material, Table S5). For the PBMC dataset with both MS and NINDs (Pooled Cohort) we selected batch of RNA-Seq library preparation and Disease-type (NINDs and MS with CIS) as covariates in the regression model. In the CD4⁺ and CD8⁺ datasets, we regressed out gender and disease type (MS and healthy) to compute the residuals (batch of RNA-Seq library preparation was corrected for in a previous step using the ComBat package) (59). For the datasets obtained from other studies of naïve B-cell, monocytes and stimulated monocytes and LCLs (Table 3), we used gender as a covariate in the regression model to compute the residuals. For the Fairfax dataset, residuals were calculated from log-transformed microarray gene expression.

Conditional analysis on LD SNPs

In the PBMC dataset for every significant SNP–gene pair association, we identified the SNPs in LD (r² >0.5) and computed the cis-eQTL beta-value and P-value for the association as described earlier in the eQTL methods. For each new significant eQTL SNP (eQTL FDR < 0.01), a stepwise association model was applied where we regressed out the effect of most-significant LD SNPs along with other covariates (batch and/or disease-type depending availability) from the expression levels (CQN values). The eQTL analysis was re-run on the remaining significant LD SNPs using the resulting residuals. In this stepwise association model, only the LD SNPs with a beta-value greater than the original MS SNP were included. To further characterize the cis-regulatory signals, we tested significant LD SNPs for the remaining eQTL effect after regressing it out with MS-associated SNP.

eQTL SNPs correlations to PBMC cell-type profile

The cell-type deconvolution method CIBERSORT (15) was applied to predict the proportion of different cell types in PBMC expression data from MS patients (Cohort 1). To understand the effect of different RNA-Seq normalization methods in predicting cell-type proportions, the CIBERSORT algorithm was first tested in RNA-Seq data obtained from CD4⁺ and CD8⁺ cells sorted from PBMCs, using the LM22 signature geneset (1000 iterations).

Subsequently, we applied the CIBERSORT to CQN-transformed expression values in the PBMC samples (Cohort 1), with the similar setting applied in sorted cell expression data. We selected only the cell types with a proportion >5% of the total PBMC mixture and aggregated them into five possible categories: CD4+ cell type (‘T cells CD4+ naïve’, ‘T cells CD4+ memory activated’ and ‘T cells CD4+ memory resting’), CD8+ cell type (‘T cells CD8+’), monocytes, B cells (‘B cells naïve’ and ‘B cells memory’) and NK cells (‘NK cells resting’ and ‘NK cells activated’). Aggregated cell-type data was corrected for batch effects using a linear model and the residual values obtained after correction were used to correlate to the 31 non-HLA SNPs and 3 HLA variants associated with significant eQTL effects in the analysis of MS PBMCs. We considered significant cell type-SNP association those with P-value <0.05 for Spearman correlation.

For SNPs that were significantly correlated to one cell type, an association model was applied, where we regressed out the effect of the corresponding cell type along with the RNA-Seq batch covariate from the expression levels (CQN values). In case of SNPs that were associated to more than one cell type a stepwise association model was implemented where corresponding cell-type proportions were added as covariates in the order of cell type-SNP association strength. The eQTL analysis was re-ran on these selected SNPs using the resulting residuals to test the independence of the SNP–gene associations from the cell type-specific effects. A potential limitation of this analysis is that the power of the dataset is not large enough to correct for any non-MS associated variants in those regions that are correlated with cell proportions.

Colocalization analysis

Bayesian-based statistical method for colocalization analysis implemented in R, ‘coloc’ package was used to test the hypothesis that independent casual variants does not exist and eQTL genes can be attributed to the MS genetic association. The method calculates the posterior probabilities for shared causal signal (H4) and independent signals (H3) from both eQTL association and disease association in a given locus (60). The test was done in loci where a minimum of 10 independent SNPs were available after LD pruning (LD r² < 0.8), within ±400-kb windows of MS-associated SNP and the prior probabilities (p1/p2 and p12) were set to the default values (1×10⁻⁴ and 1×10⁻⁵, respectively). The MS disease association P-values were obtained from Immunochip study (3) and spearman correlated P-values for eQTL genes were obtained for all the SNPs in the locus using the method described under eQTL analysis of MS-associated SNPs. For each SNP–gene the approximate bayes factor (abf) was calculated for H3 and H4 using ‘coloc.abf’ function. Then SNP–gene pairs with H3.abf greater than H4.abf were considerd as not colocalized (Supplementary Material, Table S10).

Differential expression analysis

We conducted differential gene expression analysis using DESeq2 tool (61) implemented in R. Batch effect from the RNA-Seq library preparation is added as a covariate in the design (model) and an adjusted P-value cutoff of 0.05 is used to identify significant differentially expressed genes (Supplementary Material, Table S14).

Comparison of eQTL effect sizes

To study the disease specificity, we calculated the bootstrapping derived mean of rho values of the significant eQTLs identified in Cohort 1. Then we compared estimated rho from Cohort 1 and Pooled Cohort by computing the relative change (Tables 1 and 2 and Supplementary Material, Fig. S1). For comparisons, we considered only the eQTLs with sample distribution of minimum five samples in at least two genotype-distribution (Supplementary Material, Tables S21 and S22). Calculations for Supplementary Material, Figure S1 can be found in Supplementary Material, Tables S17 and S18. To test the change in direction of the eQTL effect, we computed the mean rho values (bootstrapping derived) in a small test array dataset (Table 3) from the study (18) and relative change (ratios) were estimated from MS cohort and pooled cohort (23 MS and 7 NIND). Increase or decrease of disease enrichment (direction of change) observed in the PBMC RNA-Seq cohort was compared with the one obtained in this small test cohort and the significance of deviations from original observation (PBMC RNA-Seq cohort) was reported using binomial exact test.

To study the modification of the eQTL effect by context in sorted cells, we compared the effect ratios between stimulated and unstimulated monocytes (Fig. 5). Only effects that were significant in the reference group (unstimulated monocytes) were included in the ratio calculations. Ratios showing an effect difference of less than ±5% were depicted with white color in the corresponding heatmaps, while larger effect differences were depicted according to the color scale. Calculations for Figure 5 are found in Supplementary Material, Table S13.

Epigenomic mapping of MS-associated variants

Narrow contiguous regions of enrichment (peaks) for DNAse I hypersensitivity and Enhancer specific marks (H3K27ac and H3k4me1) based on ChromImpute (62) were obtained for PBMC and PBMC-derived immune cells (CD4⁺, CD8⁺, B cells and monocytes) from NIH Roadmap Epigenomics Mapping Consortium database (http://egg2.wustl.edu/roadmap/web_portal/index.html). Transcription factors that potentially bind to SNPs were obtained from the RegulomeDB database (63) and all these epigenetic marks were mapped to the MS-associated region.

Multiple disease association

ImmunoBase web-based tool (www.immunobase.org) was used to map other autoimmune diseases to the MS-associated SNPs.

Supplementary Material

Supplementary Material is available at HMG online.

Acknowledgements

We thank all patients who have been willing to contribute with their blood samples and make this study possible. We acknowledge International Multiple Sclerosis Genetics consortium (IMSGC) for supplying genotypes from ImmunoChip and MS replication chip used in this study. We acknowledge the Science for Life Laboratory in Stockholm for bioinformatic support in the form of a collaborative project. The computations were performed on resources provided by SNIC through Uppsala Multidisciplinary Center for Advanced Computational Science (UPPMAX) under project b2011139. We thank the Julian Knight group for sharing genotype and expression data from primary B cells and monocytes. We are grateful to Shahin Aeinehband for performing a part of the cell sorting, and to Izaura Lima Bomfim and Jenny Link for HLA imputation.

Conflict of Interest statement. T.O. has received unrestricted MS research grant, and/or honoraria for lectures and advisory boards from Biogen, Genzyme, Novartis, Merck, TEVA and Roche. All other authors declare that they have no competing interests.

Funding

This work was supported by grants from the Knut and Alice Wallenberg Foundation and the Swedish Association of Persons with Neurological Disabilities (Neuroförbundet), and Astra Zenica (AstraZeneca-Science for Life Laboratory collaboration) grant.

Author notes

The authors wish it to be known that, in their opinion, the first two authors should be regarded as joint First Authors.

The authors wish it to be known that, in their opinion, the last three authors should be regarded as joint Last Authors.