
==== Front
Proc Natl Acad Sci U S A
Proc Natl Acad Sci U S A
PNAS
Proceedings of the National Academy of Sciences of the United States of America
0027-8424
1091-6490
National Academy of Sciences

38346204
202302259
10.1073/pnas.2302259121
research-articleResearch ArticlegeneticsGenetics419
Biological Sciences
Genetics
Evidence supports a causal association between allele-specific vitamin D receptor binding and multiple sclerosis among Europeans
Adams Cameron a 1 https://orcid.org/0000-0003-2525-523X

Manouchehrinia Ali b c 1 https://orcid.org/0000-0003-4857-5762

Quach Hong L. a
Quach Diana L. a
Olsson Tomas b c d https://orcid.org/0000-0002-2938-1877

Kockum Ingrid b c d https://orcid.org/0000-0002-0867-4726

Schaefer Catherine e https://orcid.org/0000-0003-2644-055X

Ponting Chris P. f https://orcid.org/0000-0003-0202-7816

Alfredsson Lars b g h
Barcellos Lisa F. lbarcellos@berkeley.edu
a e 2
aGenetic Epidemiology and Genomics Laboratory, School of Public Health, University of California, Berkeley, CA 94720
bDivision of Neuro, Department of Clinical Neuroscience, Karolinska Institutet, Stockholm SE-171 77, Sweden
cThe Karolinska Neuroimmunology & Multiple Sclerosis Centre, Centrum for Molecular Medicine, Karolinska University Hospital, Stockholm SE-171 77, Sweden
dAcademic Specialist Center, Stockholm 113 65, Sweden
eKaiser Permanente Division of Research, Kaiser Permanente Northern California, Oakland, CA 94612
fMedical Research Council Human Genetics Unit, The Institute of Genetics and Cancer, University of Edinburgh, Western General Hospital, Edinburgh EH4 2XU, United Kingdom
gCentre for Occupational and Environmental Medicine, Region Stockholm, Stockholm 113 65, Sweden
hInstitute of Environmental Medicine, Karolinska Institutet, Stockholm SE-171 77, Sweden
2To whom correspondence may be addressed. Email: lbarcellos@berkeley.edu.
Edited by Lawrence Steinman, Stanford University, Stanford, CA; received February 26, 2023; accepted December 11, 2023

1C.A. and A.M. contributed equally to this work.

12 2 2024
20 2 2024
12 8 2024
121 8 e230225912126 2 2023
11 12 2023
Copyright © 2024 the Author(s). Published by PNAS.
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ This article is distributed under Creative Commons Attribution-NonCommercial-NoDerivatives License 4.0 (CC BY-NC-ND).

Significance

There is evidence for a causal association between 25-hydroxyvitamin D (25(OH)D) serum levels and multiple sclerosis (MS), but the underlying mechanisms are unknown. We conducted a Mendelian Randomization study to estimate causal associations between vitamin D receptor (VDR) binding at a locus and MS susceptibility using genetic instrumental variables constructed from variants previously associated with allele-specific VDR-binding (VDR-binding variant; VDR-BV) and serum 25(OH)D. Two VDR-BV (rs2881514; rs2531804) were independently associated with MS and rs2881514 demonstrated significant interaction with 25(OH)D, providing evidence for a causal association between rs2881514 and MS. Our results demonstrate that genetic variation in individual VDR binding sites contributes to MS and are relevant to other autoimmune diseases for which a role for vitamin D has been implicated.

Although evidence exists for a causal association between 25-hydroxyvitamin D (25(OH)D) serum levels, and multiple sclerosis (MS), the role of variation in vitamin D receptor (VDR) binding in MS is unknown. Here, we leveraged previously identified variants associated with allele imbalance in VDR binding (VDR-binding variant; VDR-BV) in ChIP-exo data from calcitriol-stimulated lymphoblastoid cell lines and 25(OH)D serum levels from genome-wide association studies to construct genetic instrumental variables (GIVs). GIVs are composed of one or more genetic variants that serve as proxies for exposures of interest. Here, GIVs for both VDR-BVs and 25(OH)D were used in a two-sample Mendelian Randomization study to investigate the relationship between VDR binding at a locus, 25(OH)D serum levels, and MS risk. Data for 13,598 MS cases and 38,887 controls of European ancestry from Kaiser Permanente Northern California, Swedish MS studies, and the UK Biobank were included. We estimated the association between each VDR-BV GIV and MS. Significant interaction between a VDR-BV GIV and a GIV for serum 25OH(D) was evidence for a causal association between VDR-BVs and MS unbiased by pleiotropy. We observed evidence for associations between two VDR-BVs (rs2881514, rs2531804) and MS after correction for multiple tests. There was evidence of interaction between rs2881514 and a 25(OH)D GIV, providing evidence of a causal association between rs2881514 and MS. This study is the first to demonstrate evidence that variation in VDR binding at a locus contributes to MS risk. Our results are relevant to other autoimmune diseases in which vitamin D plays a role.

multiple sclerosis
genetics
vitamin D
HHS | National Institutes of Health (NIH) 100000002 R01NS049510 Catherine SchaeferLisa Barcellos HHS | National Institutes of Health (NIH) 100000002 R01ES017080 Catherine SchaeferLisa Barcellos HHS | National Institutes of Health (NIH) 100000002 R01AI076544 Catherine SchaeferLisa Barcellos National Multiple Sclerosis Society (NMSS) 100000890 RG-19707-34664 Lisa Barcellos HHS | National Institutes of Health (NIH) 100000002 RC2AG036607 Catherine SchaeferLisa Barcellos UKRI | Medical Research Council (MRC) 501100000265 MC_UU_00007/15 Chris Ponting Magareta af Ugglas foundation NA Tomas OlssonIngrid KockumCatherine Schaefer MultipleMS project EU Horizon 2020 733161 Tomas OlssonIngrid KockumCatherine Schaefer Robert Wood Johnson Foundation (RWJF) 100000867 NA Catherine Schaefer Wayne and Gladys Valley Foundation (WGVF) 100001370 NA Catherine Schaefer The Ellison Medical Foundation NA Tomas OlssonIngrid KockumCatherine Schaefer
==== Body
pmcMultiple sclerosis (MS) is more prevalent among people residing in latitudes farther from the equator (1). This has led to the hypothesis that low vitamin D levels are associated with increased MS susceptibility (2). In observational studies, low levels of serum 25-hydroxyvitamin D (25(OH)D) (3) and lack of sunlight exposure (4, 5) are associated with an increased risk of MS, and vitamin D supplements and diets rich in vitamin D (6, 7) are associated with decreased risk of MS. However, the potential for reverse causation and unmeasured confounding factors make it difficult to assess alternative explanations. Studies employing mendelian randomization (MR), a method that utilizes genetic variation in genes associated with an exposure or other variables of interest (genetic instrumental variables or GIVs) to estimate a causal association between that exposure and disease(8), strongly support a causal relationship between low serum 25(OH)D levels and increased MS susceptibility (9–11). Many genetic variants used in these MR studies are involved in vitamin D biosynthesis, suggesting that different aspects of the vitamin D pathway, specifically transcription and expression mediated by vitamin D receptor (VDR) DNA binding, contribute to the relationship between vitamin D and MS (12, 13).

25(OH)D signals through the nuclear VDR, a ligand-regulated transcription factor that mediates all genomic actions of 25(OH)D (Fig. 1A) (14). Upon activation by vitamin D, the VDR forms the RXR/VDR heterodimer that acts as a transcription factor by binding to specific DNA regions across the genome. Most of these regions are vitamin D response elements (VDREs), a particular DNA sequence located in gene promoter regions defined by its high VDR binding affinity; however, not all VDR binding sites contain VDREs (15, 16). Through this pathway, VDR binding modulates the transcription of genes across the genome and 25(OH)D signaling through the VDR can lead to target gene expression (17). Through VDR-mediated gene transcription, vitamin D regulates calcium metabolism, cellular growth, proliferation, apoptosis, and inflammation (18). VDR binds preferentially to one allele over the other at heterozygous sites in Lymphoblastoid cell lines (LCLs) (Fig. 1B), as measured using chromatin immunoprecipitation followed by exonuclease digestion (ChIP-Exo) (19). In addition, 1000 Genomes data reveal that VDR-BVs are enriched within genomic regions associated with autoimmune conditions, including MS (19).

Fig. 1. (A) Representation of the vitamin D pathway. (B) Example of a VDR binding variant (VDR-BV) location within a VDR binding region. VDR-BVs are single-nucleotide polymorphisms associated with increased or decreased VDR binding at a VDR binding region. (C) Directed acyclic graph depicting our MR analyses. Genetic variants associated with altered VDR binding at a locus (GIVVDR-BV) are used to estimate association between VDR-Binding at a locus and MS susceptibility within our MS case–control datasets. Genetic variants from genome-wide association studies on serum vitamin D levels (GIV25(OH)D) are used to estimate the association between VDR binding and MS unbiased by horizontal pleiotropy (dashed-line).

Individual VDR-BVs have not been evaluated for association with MS susceptibility in large studies. Here, we leveraged available data on VDR-BVs and results from genome-wide association studies of 25(OH)D serum levels to construct GIVs for both VDR-BV and 25(OH)D serum levels using imputed whole-genome SNP data. GIVs can serve as proxies for exposures and other variables of interest (here for VDR-BVs and 25(OH)D) and were used in a two-sample MR study to investigate the causal relationship between VDR binding at a locus, vitamin D serum levels, and MS risk. We also utilized the dependency between VDR binding and bioavailability of 25(OH)D to develop an MR analysis framework to account for horizontal pleiotropy, a source of bias where a GIV (VDR-BV) affects disease risk outside of its effect on the variable of interest. For a true causal association, the bioavailability of 25(OH)D should modulate the effect of genetic variation in VDR binding on MS susceptibility. A significant interaction between the bioavailability of 25(OH)D and VDR binding, captured by the GIVs, would provide evidence of an association between VDR binding at a locus and MS that is not biased by horizontal pleiotropy (Fig. 1C). We hypothesized that altered VDR binding modulates transcription and expression of a target gene, thereby increasing or decreasing the risk of MS. We utilized data on individuals from three MS case–control studies and the United Kingdom Biobank. Identification of VDR-BVs associated with MS could improve the understanding of the biological mechanisms through which vitamin D acts to affect MS and further elucidate the molecular and cellular causes of MS.

Results

Characteristics of MS Cases and Controls.

The largest source of MS cases was the Swedish Human OMNI (OMNI; 6,709 MS cases and 5,881 controls), followed by the Swedish GSA (GSA; 3,718 MS cases, 1,180 controls), the UK Biobank study (UKB; 2,087 MS cases, 20,870 matched controls), and Kaiser Permanente Northern California (KPNC; 1,082 cases, 10,956 controls) (Table 1). The individuals in Swedish cohorts were genotyped on two different platforms and analyzed separately. Consistent with the female predominance in MS, all studies were primarily composed of female participants with KPNC having the highest proportion of females (80.2%) and UKB having the lowest (72.4%). The average birth year was earlier within KPNC and UKB studies (the mid-1950 s) than the other studies (1960 to 1970 s). The proportion of MS cases with at least one HLA-DRB1*15:01 allele was similar across all studies (~49 to 57%), with UKB cases having the lowest carriage rate. All analyses were performed separately within KPNC, GSA, OMNI, and UKB and combined with meta-analysis.

Table 1. Characteristics of MS cases and controls

	KPNC	UKB	GSA	OMNI	
	Case	Control	Case	Control	Case	Control	Case	Control	
n	1,084	10,956	2,087	20,870	3,718	1,180	6,709	5,881	
Sex, n (%)									
  Male	215 (19.8)	2,165 (19.8)	576 (27.6)	5,760 (27.6)	1,031 (27.7)	342 (29.0)	1,808 (26.9)	1,455 (24.7)	
  Female	869 (80.2)	8,791 (80.2)	1,511 (72.4)	15,110 (72.4)	2,687 (72.3)	838 (71.0)	4,901 (73.1)	4,426 (75.3)	
Birth year, mean (SD)	1957.6 (8.9)	1957.1 (8.6)	1952.7 (7.6)	1952.7 (7.6)	1973.1 (12.4)	1968.9 (13.7)	1961.0 (13.7)	1960.5 (13.3)	
DRB1*15:01 carrier, n (%)									
  0 alleles	506 (46.7)	8,280 (75.6)	1,050 (50.3)	15,271 (73.2)	1,621 (43.6)	843 (71.4)	2,913 (43.4)	4,052 (68.9)	
  1 or 2 alleles	578 (53.3)	2,676 (24.4)	1,037 (49.7)	5,599 (26.8)	2,097 (56.4)	336 (28.5)	3,796 (56.6)	1,829 (31.1)	
25(OH)D GIV, Jiang et al.* [mean (SD)]	0.285 (0.07)	0.288 (0.07)	0.292 (0.07)	0.293 (0.07)	0.289 (0.07)	0.293 (0.07)	0.286 (0.07)	0.290 (0.07)	
25(OH)D GIV, Revez et al.* [mean (SD)]	2.150 (0.12)	2.163 (0.19)	2.184 (0.19)	2.187 (0.19)	2.057 (0.19)	2.074 (0.18)	2.066 (0.19)	2.070 (0.19)	
Abbreviations: GSA, Swedish cohort genotyped on global screening array; Human Omni, Swedish cohort genotyped on human omni express array; KPNC, Kaiser Permanente Northern California MS case–control study; GIV, Genetic Instrumental Variable; UKB, UK Biobank.

*GIVs calculated from independent GWAS variants (Linkage disequilibrium R2 < 0.001) with P < 5 × 10−8.

Linkage Disequilibrium between VDR-BVs and MS GWAS Risk Variants.

Within the European samples from the 1000 Genomes project, linkage disequilibrium (LD) between VDR-BVs and established MS GWAS variants was generally low. Of the 112 VDR-BVs considered in this investigation, only two VDR-BVs had an r2 > 0.2 with a known MS GWAS variant (VDR-BV: rs55792977, MS GWAS: rs13385171, r2 = 0.48; VDR-BV: rs13098781, MS GWAS: rs9863496, r2 = 0.34), indicating no correlation between MS GWAS risk variants and nearly all VDR-BV variants.

25(OH)D Serum Levels Are Casually Associated with MS Susceptibility.

Two GIVs constructed using summary statistics from two recent GWASs on 25(OH)D serum levels by Jiang et al. and Revez et al. (hereafter, GIV25OHD; see Materials and Methods) were used in our analyses (12, 13). The mean GIV25(OH)D for increased 25(OH)D serum levels was lower among MS cases compared to controls in all studies (Table 1). The difference for both the Jiang et al. GIV25(OH)D and Revez et al. GIV25(OH)D means between MS cases and controls was largest among Sweden participants genotyped using the GSA chip. The difference in the means was smallest across both GIV25(OH)D between UKB MS cases and controls. Consistent with previous research, a decreased 25(OH)D serum level captured by each GIV25(OH)D was associated with an increased risk of MS (Fig. 2 A and B). The magnitude of the association was greater for the Jiang et al. GIV25(OH)D (OR: 1.85, 95% CI: 1.30–2.63) compared to the Revez et al. GIV25(OH)D (OR: 1.32, 95% CI: 1.07–1.61), likely owing to the difference in the number of SNPs within each GIV (Materials and Methods) and reported effect sizes. Estimates for between study variance were low for the Jiang et al. GIV25(OH)D (I2 = 0%, PQ = 0.56) indicating associations were consistent across all studies. While there was some evidence for moderate heterogeneity for the Revez et al. GIV25(OH)D (I2 = 55%, PQ = 0.09), the direction of effect was consistent across all studies, providing support for a uniform pattern across the studies, and strong evidence for causal association between low 25(OH)D serum level and MS risk.

Fig. 2. Results from meta-analyses of association between 25(OH)D genetic instrumental variables (GIV25OHD) and VDR-BV rs2881514 and MS. GIVs25OHD calculated using summary statistics from genome-wide association studies on serum 25(OH)D levels. (A) GIVs25OHD calculated using summary statistics from Jiang et al. (13); (B) GIVs25OHD calculated using summary statistics from Revez et al. (12). (C) Meta-analysis results for association between GIV25(OH)D rs2881514 and MS and (D) evidence of interaction between rs2881514 and Jiang et al. 25(OH)D GIV. τ2 = estimate of between study heterogeneity; I2 = proportion of variance due to heterogeneity; and P = P-value for Cochran’s Q test of heterogeneity.

VDR-BV Associations with MS.

GIVs for each VDR-BV (hereafter, GIVVDR) were constructed (Materials and Methods). Of the 112 GIVsVDR included in our analyses, seven GIVsVDR were associated with MS at P < 0.05 (Table 2). Following correction for multiple testing, two GIVVDR remained significant (PFDR < 0.05). Meta-analyses of these two GIVsVDR demonstrated evidence of an association between variation within VDR binding and MS susceptibility across all four studies. GIVVDR rs2881514 was associated with an increased risk of MS [Meta-analysis Odds Ratio (OR): 1.10, 95% CI: 1.05–1.15, P = 9.4×10-5 ; Fig. 2C] and GIVVDR rs2531804 was associated with a decreased risk of MS (Meta-analysis OR: 0.82, 95% CI: 0.73–0.92, P = 6.4×10-4 ). There was evidence of interaction between 13 GIVsVDR and GIV25OHD at P < 0.05 (SI Appendix, Table S1). No interactions were significant after correcting for multiple testing; however, GIVVDR rs2881514, the top finding from the independent models, was among these 13 (GIVVDR rs2881514 × Jiang et al. GIV25OHD Meta-analysis OR: 2.17, 95% CI: 1.10–4.29, P = 0.025; Fig. 2D). The direction of this association provided evidence of interaction between increased VDR binding at GIVVDR rs2881514 and decreased bioavailability of 25(OH)D measured with Jiang et al. GIV25OHD conferring increased risk of MS.

Table 2. VDR binding variants (VDR-BVs) instrumental variable (IV) associated with MS susceptibility at P < 0.05

GIVVDR*	chr:bp	GREAT annotation†	OR‡	95% CI	P	P FDR ‡	
rs2881514_A	chr3:16553679	RFTN1 (+1533); OXNAD1 (+246974)	1.10	1.05–1.15	9.4E-05	0.011	
rs2531804_A	chr6:28411302	ZSCAN23 (-59)	0.82	0.73–0.92	6.4E-04	0.036	
rs12048389_T	chr1:107538722	PRMT6 (-60578)	0.92	0.87–0.99	1.8E-02	0.660	
rs55792977_T	chr2:65650863	SPRED2 (+8447); ACTR2 (+195893)	1.05	1.00–1.09	3.2E-02	0.660	
rs7309003_C	chr12:97751789	NEDD1 (+450546)	0.95	0.91–1.00	3.2E-02	0.660	
rs10995246_C	chr10:64391845	ADO (-172670); ZNF365 (+257895)	1.06	1.00–1.13	3.5E-02	0.660	
rs10232857_C	chr7:55601335	VOPP1 (+38882); LANCL2 (+168195)	0.93	0.87–1.00	4.9E-02	0.780	
Abbreviations: bp, base pair; chr; chromosome; CI, Confidence Interval; FDR, false discovery rate; GIV, Genetic Instrumental Variable; GREAT, Genomic Regions Enrichment of Annotations Tool; OR, Odds Ratio; VDR-BV, Vitamin D Receptor Binding Variant.

*VDR-BV using as instrumental variable for allelic specific binding. Allele included is the allele associated with increased binding.

†Distance (bp) from VDR-BV to TSS of the nearest upstream and downstream gene from rGREAT.

‡Odds ratios are from random effect meta-analyses combining estimates from the four studies. ORs for association between VDR-BV and MS.

VDR-BVs with Evidence of MS Association Are eQTLs and sQTLs.

Several of the GIVVDR that were independently associated at P < 0.05 with MS or showed evidence of interaction with GIV25(OH)D are expression quantitative trait loci (eQTL) or splicing quantitative trait loci (sQTL) for one or more human tissues in GTEx (SI Appendix, Fig. S1) (20). GIVVDR rs2881514 is an eQTL for RFTN1 in esophageal tissue and an sQTL for RFTN1 in EBV-transformed lymphocytes and thyroid tissues. RFTN1 encodes Raftlin, a B cell-specific major raft protein necessary for lipid raft integrity and B-cell receptor signal transduction (21). Other GIVVDR, including rs2531804, rs55811049, rs961320, rs871699, and rs62200158, which showed evidence of association with MS at P < 0.05, are eQTLs or sQTLs in brain or skin tissue, among others.

Discussion

This is the first study to observe evidence for associations between variation in VDR binding at a locus and MS susceptibility, providing further evidence of the important role of the vitamin D pathway in MS. Our results support the hypothesis that alterations in VDR binding and subsequent modulations in vitamin D gene transcription and expression contributes to developing MS. We found two GIVVDR, rs2881514 and rs2531804, were independently associated with MS susceptibility. There was also evidence of interaction between several GIVsVDR and GIV25(OH)D and MS risk. Although no observed interactions were significant after correction for multiple testing, there was evidence of interaction between GIVVDR rs2881514 and the Jiang et al. GIV25(OH)D at P < 0.05. It is important to note that detecting interactions requires substantially more statistical power than main effects. Regardless, evidence of an association between a GIVVDR and MS in independent and interaction models is promising (22). Both GIV25(OH)D for 25(OH)D serum levels were significantly associated with MS; specifically, lower 25(OH)D serum was associated with increased risk of MS in all four studies.

There are causal associations between vitamin D serum levels and MS risk using MR methods (9–11). Our findings using GIVs constructed using summary statistics from two recent GWASs on 25(OH)D serum levels are consistent with previous MR studies and further confirm vitamin D insufficiency as an important risk factor for MS. Measures of heterogeneity indicate little variation between studies for Jiang et al. GIV25(OH)D, although there was some evidence of moderate heterogeneity for the Revez et al. GIV25(OH)D (I2 = 55%, PQ = 0.09), which might be explained by the larger number of SNPs included in the Revez et al. GIV25(OH)D. Associations were modest for both GIV25(OH)D among UKB participants. Unlike the other studies, MS cases from UKB did not have their MS status confirmed by a neurologist. Instead, diagnoses were inferred from electronic health records or self-report and may have resulted in a small amount of misclassification of MS case status in this cohort. Evidence of this may be seen in the lower frequency of carriage of HLA-DRB1*15:01, the major risk allele for MS risk, in UKB MS cases (49.7%) compared to cases in the other three studies (>53%).

The top GIVVDR associated with MS from meta-analyses was rs2881514. The A allele of this SNP was associated with increased VDR binding in LCLs. Evidence from results within each study and meta-analyses indicate a harmful effect of increased VDR binding at that locus on MS risk. Further, the interaction observed between GIVVDR rs2881514 and GIV25(OH)D demonstrates evidence for interaction between increased VDR binding at rs2881514 and decreased 25(OH)D serum levels, yielding a substantial increase in the risk of MS (Meta-analysis OR: 2.17, 95% CI: 1.10–4.29). This association was consistent among UKB, OMNI, and GSA MS cases and controls but was slightly protective in KPNC, though very wide confidence intervals were observed (KPNC OR: 0.79, 95% CI: 0.14–4.47). GIVVDR rs2881514 is located on chromosome 3, 1,533 bases downstream of the transcription start site of RFTN1. In GTEx, rs2881514 is an eQTL for RFTN1 in esophageal tissue and an sQTL for RFTN1 in EBV-transformed lymphocytes and thyroid tissues, providing evidence for a role of rs2881514 in RFTN1 expression. RFTN1 encodes Raftlin, a protein that is critical for producing lipid rafts, which are membrane microdomains that play a crucial role in B cell activation through B cell receptor signaling (21). Raftlin is essential for clathrin-dependent endocytosis of Toll-like Receptor (TLR) 3 ligand in human epithelial cells and myeloid dendritic cells (23), and for LPS-induced TLR4 internalization and Toll-IL-1R domain-containing adaptor molecule-1 (TICAM-1) signaling in human monocyte-derived DCs and macrophages (24). TLR3 and TLR4 signaling pathways have been implicated in MS pathogenesis and symptom modulation (25). In mice, experimental autoimmune encephalomyelitis (EAE), a Th1-/Th17-mediated disease, was more severe among mice with wild-type Raftlin expression compared to Raftlin-deficient mice (23). RFTN1 has not been directly implicated in MS risk; however, there is increased expression of RFTN1 in chronic active brain lesions from MS cases compared to healthy brain tissue from controls (26). More work is needed to understand the relationship between RFTN1 and vitamin D, and its role in MS.

Several of the GIVVDR associated with MS had evidence of SNP-associated tissue-specific expression and alternative splicing in GTEx (SI Appendix, Fig. S1). rs25318104, rs9621320, rs871699, rs689384, rs62200158, rs558110449, and rs2286576 show evidence of being an eQTL or sQTL in one or more brain tissues. Many of these SNPs are eQTLs or sQTLs across various tissues. Of the set of VDR-BVs tested as GIVVDR in the current study, those with evidence of association with MS susceptibility were enriched for genes involved in immune-related processes (leukocyte-mediated immunity, immune effector process, myeloid leukocyte activation, and phospholipase activity), vesicles, and exocytosis (SI Appendix, Table S2). Annotations to eQTLs and sQTLs and the enrichment results provide evidence supporting a mediating role of allele-specific VDR binding in the expression of nearby genes and immune pathways in MS risk.

Observational studies have shown that increased dietary intake of vitamin D is associated with reduced risk of MS, but it is unknown whether interventional vitamin D supplementation would yield protective effects on disease susceptibility (6, 27). Studies investigating the effects of vitamin D supplementation on MS progression have found associations between higher serum levels of vitamin D and reduced hazard of relapse and a reduction in incidence of T2 lesions (28, 29). Our results show how genetic variation in the vitamin D pathway affects MS disease risk. Future studies investigating the protective effects of vitamin D supplementation on disease risk and progression should consider genetic variation in vitamin D pathway genes as part of their design.

This study had several strengths. It is the first to investigate the association of VDR binding at a locus with MS among cases and controls. An investigation of variation in VDR binding among cases and controls would not be feasible in large sample sizes without utilizing MR methodology. We used previously identified VDR-BVs to create GIVs for VDR binding at a locus (19). Another strength was the large sample size and use of data from independent studies. There were also some limitations, including the possibility of horizontal pleiotropy, biasing associations between GIVVDR and MS. Each GIVVDR was constructed using only one VDR-BV, preventing the use of standard sensitivity analyses for assessing the validity of MR assumptions in our VDR-BV analyses. However, our use of GIVs for 25(OH)D serum level in the analyses allowed us to estimate the association between GIVVDR and MS unbiased by horizontal pleiotropy. Additionally, all but two of the VDR-BV were independent of known MS GWAS variant risks. MS cases and controls in the study were of European ancestry to reduce the impact of population stratification; however, this approach limits the generalizability of our findings to non-European populations. Further, VDR-BVs were identified in LCLs which are EBV-transformed B-cell lines. Variation in VDR binding associated with MS likely occurs in other lymphocytes, including CD4+ and CD8+ T cells and it is unknown whether the VDR-BVs active in LCLs are also active in these T cells. We were also only able to assess variation in VDR binding at 112 loci, which represent a small portion of genomic regions where VDR binds. Lastly, there are likely to be mechanisms other than VDR binding that underlie the observed VDR-MS associations reported here, including alternative splicing, effect of microRNAs, and combinatorial protein binding. Our study was designed to identify causal mechanisms of MS susceptibility that are mediated by altered VDR binding due to DNA variation.

This study demonstrates that genetic variation in VDR binding at a single locus contributes to MS susceptibility. Our results highlight the importance of the vitamin D pathway in MS pathogenesis. Our results are also relevant to other autoimmune and inflammatory diseases and several cancers for which a role for vitamin D has been suggested. Future studies of VDR binding and MS should identify VDR-BVs in lymphocytes not captured by LCLs, including CD4+ and CD8+ T cells.

Materials and Methods

MS Case–Control Studies.

Individual-level case–control data used in this research was from one USA-based case–control study and two Sweden-based MS case–control studies. MS cases and controls from the US were from the Kaiser Permanente Northern California (KPNC) MS Research Program (30, 44). Additional KPNC controls were participants of the Genetic Epidemiology Research on Adult Health and Aging (GERA; dbGaP phs000674.v2. p2) (30). From Sweden, MS cases and controls were from the Nationwide Epidemiological Investigation of Multiple Sclerosis (EIMS), the Genes and Environment in Multiple Sclerosis (GEMS), IMSE study of the effect of immunomodulatory drugs, and the Karolinska Hospital STOPMS study (32, 33, 45). Methods for identification and confirmation of MS diagnosis in each study are provided within the corresponding publications. Briefly, MS cases from KPNC, EIMS, GEMS, IMSE, and STOPMS had their disease status confirmed by an independent neurologist. All study participants provided written informed consent, and all studies obtained approval from the Institutional Review Boards of KPNC, local Ethical Committees, and the University of California, Berkeley.

UK Biobank.

The UK Biobank (UKB, http://www.ukbiobank.ac.uk) is a prospective cohort study of approximately 500,000 individuals from the United Kingdom (31). Recruitment took place between 2006 and 2010 in 22 assessment centers located across the United Kingdom. Participants were aged 40 to 69 at the time of recruitment. MS cases were identified using UKB field 31043-0.0 “Source of report of G35 (MS).” Controls were defined as those not reporting MS symptoms or other demyelinating diseases at the time of enrollment and were frequency matched to cases by year of birth ( ± 2 y) and sex at a ratio of 10 controls per MS case. Participants included in analyses were unrelated individuals of European ancestry with imputed genetic data and no sex chromosome aneuploidies. Those who had withdrawn consent or were recommended for genetic analysis exclusion were excluded. Data were accessed under the approval of UKB within project 69668 (46). All participants gave prior written informed consent, and the study was conducted following the principles of the Declaration of Helsinki.

Genotype and Exposure Assessment.

All participants in the MS case–control studies (KPNC, GEMS, EIMS, IMSE, and STOP-MS) completed an interview or self-reported questionnaire related to MS disease events, reproductive history, and environmental exposures (30, 32–34). Participants provided blood or saliva samples for genotyping. SNP genotyping was performed using the Illumina Infinium 660 K BeadChip Array and Human Omni Express Array (KPNC), Illumina Global Screening Array or Human Omni Expression (GEMS, EIMS, IMSE, and STOPMS), and Affymetrix Axiom Array (GERA). The data from GEMS, EIMS, IMSE, and STOPMS were analyzed in two cohorts depending on which genotyping array was used: Human Omni Express (OMNI) or Global screening array (GSA). Details of genotyping and imputation for KPNC (9), GSA AND OMNI (35), and UKB have been described previously. Participants with missing genotypes that met QC thresholds (info score > 0.8, missingness per SNP < 0.05, missingness per cohort < 0.05, and minor allele frequency (MAF) > 0.05) were imputed using the mean MAF within each study.

Genetic Instruments for VDR Binding.

MR is a method that uses one or more genetic variants associated with an exposure or other variable of interest as a GIV to estimate causal associations between that exposure and an outcome. GIVs are calculated by summing the weighted contributions of individual SNPs, based on their effect sizes derived from GWASs or other genetic studies. In this study, GIVVDR-BV for each individual were constructed using the following equation:GIVVDR-BVij=dosageij ×βVDR-BVj,

where dosageij is the number of effect alleles at VDR-BVj carried by individual i and βVDR-BVi is the effect size for the association between VDR-binding at a locus and alleles at VDR-BVj . There are three primary MR assumptions: the variant(s) must be associated with the exposure, there are no unmeasured confounders of the association between the genetic variant(s) and the outcome, and the genetic variants only affect the outcome through the exposure (8). The occupancy of VDR at a given genomic locus was the exposure of interest in this study. Our objective was to use SNPs associated with VDR binding at a locus previously identified by Gallone et al (19) as GIVs to identify associations between variation in VDR binding at a locus and MS susceptibility. Occupancy of VDR was quantified genome-wide at heterozygous sites using ChIP-exo sequencing read data from sixteen genotyped calcitriol-stimulated LCLs. To satisfy the first assumption of MR, at least one variant statistically associated with altered VDR occupancy at a locus was required (8). SNPs associated with VDR binding in cis were identified by observing differential VDR occupancy over each allele of a given heterozygous SNP. Allele-specific VDR binding was identified using AlleleSeq, modeling reads mapping to each allele using the binomial distribution and applying a two-tailed test (36). VDR-BVs were called at heterozygous sites with five or more ChIP-Exo reads, with at least read being mapped to each allele in order to reduce unwanted effects of PCR duplication. After accounting for copy number variation (33) and correction for multiple testing using AlleleSeq’s computational simulation approach (false discovery rate < 2%), 305 VDR-BVs were identified. Of these, 112 VDR-BVs were present in all imputed case–control datasets for analysis. SNP-VDR binding effect sizes were estimated by regressing the read counts from ChIP-exo data at a heterozygous VDR binding region described above against the corresponding VDR-BV alleles. For each VDR-BV, the GIV was derived by multiplying the SNP-VDR binding effect size by the number alleles carried by each participant (hereafter, GIVVDR).

Pleiotropy and Bioavailability of Vitamin D.

The causal effect estimates from MR analysis are unbiased in the absence of pleiotropic effects. A common approach to estimating the pleiotropic bias is based on over-dispersion across multiple instrumental SNPs for a given exposure. However, sensitivity analyses for multi-SNP MR instruments were not applicable as there was only one instrumental SNP for each GIVVDR. Given that VDR binding is dependent upon the bioavailability of vitamin D, we used variants and summary statistics from two recent GWAS on serum 25(OH)D to construct GIVs for bioavailability of vitamin D (see below) to perform subgroup analyses (12, 13). Interaction between a GIVVDR and bioavailability of 25(OH)D was considered evidence of an association between VDR binding and MS susceptibility not biased by horizontal pleiotropy.

Genetic Instruments for 25(OH)D.

SNPs and the estimated effect sizes from two recent GWAS on serum 25(OH)D were used to construct each GIV25OHD (12, 13). GWAS summary statistics were extracted from MR-Base R platform (37). The extract_instruments() function in the TwoSampleMR R package was used to identify independent variants with GWAS P < 1 × 10−8 within genomic windows of 10,000 Kb (R2 < 0.001) using a European LD reference panel from 1000 Genomes Phase 3. These procedures yielded eight SNPs for the Jiang et al. GIV25(OH)D and 103 SNPs for the Revez et al. GIV25(OH)D (SI Appendix, Tables S3 and S4). After identifying these SNPs, each GIV25(OH)D was calculated in Plink1.9 using --score (38) which uses the following equation:GIV25OHDi=∑kKβk×dosageik,

where dosageik is the number of effect alleles at GWAS variant k carried by individual i and βk is the effect size for the association between reduced 25(OH)D serum levels at and GWAS variant k.

Statistical Methods.

Logistic regression was used to estimate the associations between all GIVs (GIV25OHD and GIVVDR) and MS susceptibility within each study:logitPMS=1∼ β0+β1GIV25OHD,

logitPMS=1∼ β0+β1GIVVDRj.

Additionally, to estimate an association not biased by horizontal pleiotropy, we estimated the interaction between each GIVVDR and GIV25OHD for MS susceptibility:logitPMS=1∼β0+β1GIVVDRj+β2GIV25OHD+β3GIVVDRj×GIV25OHD,

where β3 is the estimate of multiplicative interaction between VDR binding at a locus and the bioavailability of 25(OH)D. A significantly non-zero value of β3 indicates that pleiotropy does not explain all the observed association and therefore, the association between VDR-BV and MS is causal.

All models were adjusted for sex (male or female), quintiles of birth year, carriage of the HLA-DRB1*15:01 allele, and the first six genome-wide principal components. Random effects meta-analysis was used to combine study-specific associations. All meta-analyses were performed using the metagen() function in the Meta R package (39). The DerSimonian-Laird estimator was used to estimate the between-study variance. Between-study heterogeneity was assessed using Cochran’s Q statistic and Higgins & Thompson’s I2 statistic. P-values from associations between each GIV25OHD and MS and the interaction parameter between each GIVVDR and GIV25OHD were corrected for multiple tests using the Benjamini–Hochberg method (40).

Annotation of VDR-BVs to the nearest gene transcription start sites (TSS) and GO enrichment analysis was performed using rGREAT (41). Tissue-specific expression quantitative trait loci (eQTL) and splicing quantitative trait loci (sQTL) for VDR-BVs were from the GTEx v8 Project (20). LD between VDR-BVs and MS GWAS risk variants was estimated using the ld_matrix() function from the ieugwasr R package using the European LD reference panel from 1000 Genomes Phase 3 (42).

All statistical analyses were conducted using R 4.0.2 (43, 47).

Supplementary Material

Appendix 01 (PDF)

This work was supported in part by the following: NIH NINDS R01NS049510 (PI Barcellos), NIH NIEHS R01ES017080 (PI Barcellos), NIH NIAID R01AI076544 (PI Barcellos), National MS Society RG-19707-34664 (PI Barcellos), NIH RC2AG036607 (PI Schaefer), UK Medical Research Council MC_UU_00007/15 (PI Ponting), the Robert Wood Johnson Foundation (PI Schaefer), the Wayne and Gladys Valley Foundation (PI Schaefer), the Ellison Medical Foundation (PI Schaefer), the Magareta af Ugglas foundation (PI Olsson), and MultipleMS project EU Horizon 2020 (project 733161, PI Kockum). We also thank Chris Cotsapas and Neil Clark for their contributions to this work.

Author contributions

C.A., A.M., I.K., C.P.P., L.A., and L.F.B. designed research; C.A., A.M., H.L.Q., and D.L.Q. performed research; C.A. and A.M. analyzed data; and C.A., A.M., T.O., I.K., C.S., C.P.P., L.A., and L.F.B. wrote the paper.

Competing interests

The authors declare no competing interest.

Data, Materials, and Software Availability

Metadata for KPNC (https://doi.org/10.6084/m9.figshare.22145966.v1) (44) as well as GSA and OMNI participants (https://doi.org/10.6084/m9.figshare.22047314.v1) (45) is available. Complete KPNC data can be accessed by contacting Lisa Barcellos (lbarcellos@berkeley.edu) and Lynn Hollyer (lhollyer@berkeley.edu). Requests will be reviewed by the IRB, and data can be shared upon approval. Complete GSA and OMNI data can be accessed by contacting Ingrid Kockum (Ingrid.Kockum@ki.se). Genotype and phenotype data for UKB participants can be accessed from UKB (https://www.ukbiobank.ac.uk/enable-your-research/about-our-data/genetic-data) (46). GERA data are available on dbGaP (phs000674.v2.p2). R scripts for regression analyses and meta-analyses are available here: https://doi.org/10.5281/zenodo.7647412 (47). Summary statistics from Gallone et al. for VDR-BVs used in this study (SI Appendix, Table S5) and study-specific and meta-analysis associations between GIVVDR-BV and MS susceptibility for all VDR-BVs are available in supplement (SI Appendix, Tables S6–S9).

Supporting Information

This article is a PNAS Direct Submission.
==== Refs
1 G. Rosati, The prevalence of multiple sclerosis in the world: An update. Neurol. Sci. 22 , 117–139 (2001).11603614
2 P. Goldberg, Multiple sclerosis: Vitamin D and calcium as environmental determinants of prevalence. Int. J. Environ. Stud. 6 , 19–27 (1974).
3 K. L. Munger, L. I. Levin, B. W. Hollis, N. S. Howard, A. Ascherio, Serum 25-hydroxyvitamin D levels and risk of multiple sclerosis. JAMA 296 , 2832 (2006).17179460
4 I. A. F. van der Mei , Past exposure to sun, skin phenotype, and risk of multiple sclerosis: Case–control study. BMJ 327 , 316 (2003).12907484
5 M. T. Kampman, T. Wilsgaard, S. I. Mellgren, Outdoor activities and diet in childhood and adolescence relate to MS risk above the Arctic Circle. J. Neurol. 254 , 471–477 (2007).17377831
6 K. L. Munger , Vitamin D intake and incidence of multiple sclerosis. Neurology 62 , 60–65 (2004).14718698
7 A. Ascherio, K. L. Munger, K. C. Simon, Vitamin D and multiple sclerosis. Lancet Neurol. 9 , 599–612 (2010).20494325
8 G. D. Smith, S. Ebrahim, “Mendelian randomization”: Can genetic epidemiology contribute to understanding environmental determinants of disease? Int. J. Epidemiol. 32 , 1–22 (2003).12689998
9 B. Rhead , Mendelian randomization shows a causal effect of low vitamin D on multiple sclerosis risk. Neurol. Genet. 2 , e97 (2016).27652346
10 L. E. Mokry , Vitamin D and risk of multiple sclerosis: A mendelian randomization study. PLoS Med. 12 , e1001866 (2015).26305103
11 M. A. Gianfrancesco , Evidence for a causal relationship between low vitamin D, high BMI, and pediatric-onset MS. Neurology 88 , 1623–1629 (2017).28356466
12 J. A. Revez , Genome-wide association study identifies 143 loci associated with 25 hydroxyvitamin D concentration. Nat. Commun. 11 , 1–12 (2020).31911652
13 X. Jiang , Genome-wide association study in 79,366 European-ancestry individuals informs the genetic architecture of 25-hydroxyvitamin D levels. Nat. Commun. 9 , 260 (2018).29343764
14 D. D. Bikle, Vitamin D metabolism, mechanism of action, and clinical applications. Chem. Biol. 21 , 319–329 (2014).24529992
15 S. V. Ramagopalan , A ChIP-seq defined genome-wide map of vitamin D receptor binding: Associations with disease and evolution. Genome Res. 20 , 1352–1360 (2010).20736230
16 P. Tuoresmäki, S. Väisänen, A. Neme, S. Heikkinen, C. Carlberg, Patterns of genome-wide VDR locations. PLoS ONE 9 , e96105 (2014).24787735
17 J. W. Pike, M. B. Meyer, The vitamin D receptor: New paradigms for the regulation of gene expression by 1,25-dihydroxyvitamin D3. Endocrinol. Metab. Clin. North Am. 39 , 255–269 (2010).20511050
18 T. K. Wöbke, B. L. Sorg, D. Steinhilber, Vitamin D in inflammatory diseases. Front. Physiol. 5 , 244–244 (2014).25071589
19 G. Gallone , Identification of genetic variants affecting vitamin D receptor binding and associations with autoimmune disease. Hum. Mol. Genet. 26 , 2164–2176 (2017).28335003
20 The GTEx Consortium, The GTEx Consortium atlas of genetic regulatory effects across human tissues. Science 369 , 1318–1330 (2020).32913098
21 K. Saeki, The B cell-specific major raft protein, Raftlin, is necessary for the integrity of lipid raft and BCR signal transduction. EMBO J. 22 , 3015–3026 (2003).12805216
22 W. J. Gauderman , Update on the state of the science for analytical methods for gene–environment interactions. Am. J. Epidemiol. 186 , 762–770 (2017).28978192
23 A. Watanabe , Raftlin Is involved in the nucleocapture complex to induce poly(I:C)-mediated TLR3 activation. J. Biol. Chem. 286 , 10702–10711 (2011).21266579
24 M. Tatematsu , Raftlin controls lipopolysaccharide-induced TLR4 internalization and TICAM-1 signaling in a cell type-specific manner. J. Immunol. 196 , 3865–3876 (2016).27022195
25 C. Zheng, J. Chen, F. Chu, J. Zhu, T. Jin, Inflammatory role of TLR-MyD88 signaling in multiple sclerosis. Front. Mol. Neurosci. 12 , 314 (2020).31998072
26 D. A. E. Hendrickx , Gene expression profiling of multiple sclerosis pathology identifies early patterns of demyelination surrounding chronic active lesions. Front. Immunol. 8 , 1810 (2017).29312322
27 M. B. Sintzel, M. Rametta, A. T. Reder, Vitamin D and multiple sclerosis: A comprehensive review. Neurol. Therapy 7 , 59–85 (2018).
28 S. Simpson , Higher 25-hydroxyvitamin D is associated with lower relapse risk in multiple sclerosis. Ann. Neurol. 68 , 193–203 (2010).20695012
29 E. M. Mowry , Vitamin D status predicts new brain magnetic resonance imaging activity in multiple sclerosis. Ann. Neurol. 72 , 234–240 (2012).22926855
30 Y. Banda , Characterizing race/ethnicity and genetic ancestry for 100,000 subjects in the genetic epidemiology research on adult health and aging (GERA) cohort. Genetics 200 , 1285–1295 (2015).26092716
31 C. Sudlow , UK Biobank: An open access resource for identifying the causes of a wide range of complex diseases of middle and old age. PLoS Med. 12 , e1001779 (2015).25826379
32 F. B. S. Briggs , Adverse socioeconomic position during the life course is associated with multiple sclerosis. J. Epidemiol. Commun. Health 68 , 622–629 (2014).
33 A. K. Hedström, J. Hillert, T. Olsson, L. Alfredsson, Alcohol as a modifiable lifestyle factor affecting multiple sclerosis risk. JAMA Neurol. 71 , 300 (2014).24395432
34 A. Hedström, T. Åkerstedt, T. Olsson, L. Alfredsson, Shift work influences multiple sclerosis risk. Mult. Scler. J. 21 , 1195–1199 (2015).
35 S. Baranzini, S. Sawcer, International Multiple Sclerosis Genetics Consortium, MultipleMS Consortium, Genetic analysis of multiple sclerosis severity identifies a novel locus and implicates CNS resilience as a major determinant of outcome (2022) 10.21203/rs.3.rs-1723574/v1.
36 J. Rozowsky , AlleleSeq: Analysis of allele-specific expression and binding in a network framework. Mol. Syst. Biol. 7 , 522 (2011).21811232
37 G. Hemani , The MR-Base platform supports systematic causal inference across the human phenome. eLife 7 , e34408 (2018).29846171
38 C. C. Chang , Second-generation PLINK: Rising to the challenge of larger and richer datasets. Gigascience 4 , 7 (2015).25722852
39 S. Balduzzi, G. Rücker, G. Schwarzer, How to perform a meta-analysis with R: A practical tutorial. Evid. Based Ment. Health 22 , 153–160 (2019).31563865
40 D. Yekutieli, Y. Benjamini, The control of the false discovery rate in multiple testing under dependency. Ann. Stat. 29 , 1165–1188 (2001).
41 Z. Gu, D. Hübschmann, rGREAT: an R/bioconductor package for functional enrichment on genomic regions. Bioinformatics 39 , btac745 (2023).
42 B. Elsworth , The MRC IEU OpenGWAS data infrastructure. bioRxiv [Preprint] (2020). 10.1101/2020.08.10.244293 (Accessed 13 May 2022).
43 R Core Team, R: A language and environment for statistical computing (R Foundation for Statistical Computing, Vienna, Austria, 2021). https://www.R-project.org/. Accessed 22 February 2022.
44 C. Adams, KPNC Metadata. Figshare. 10.6084/m9.figshare.22145966.v1. Deposited 23 February 2023.
45 I. Kockum, FigShare_VitD.xlsx. Figshare. 10.6084/m9.figshare.22047314.v1. Deposited 8 February 2023.
46 C. Bycroft , The UK Biobank resource with deep phenotyping and genomic data. UK Biobank. https://www.ukbiobank.ac.uk/enable-your-research/about-our-data/genetic-data. Accessed 6 July 2021.
47 C. Adams, adams-cam/MS_VDR_MR: v1. Zenodo. 10.5281/zenodo.7647412. Deposited 16 February 2023.
