Molecular mechanism underlying the impact of vitamin D on disease activity of MS.

OBJECTIVE: Some previous studies suggest modest to strong effects of 25-hydroxyvitamin D (25(OH)D) on multiple sclerosis (MS) activity. The objective of this study was to explore the mechanistic rationale that may explain potential clinical effects of 25(OH)D.
METHODS: This study measured serum 25(OH)D levels and global gene expression profiles over a course of up to 2 years in patients starting treatment with interferon beta-1b (IFNB-1b) after a clinically isolated syndrome. MS disease activity was assessed by the number of gadolinium-enhancing lesions present on repeated magnetic resonance imaging (MRIs).
RESULTS: The number of gadolinium-enhancing lesions was highly significantly associated with 25(OH)D levels. Conducting various systems-level analyses on the molecular level, multiple lines of evidence indicated that 25(OH)D regulates expression dynamics of a large gene-gene interaction system which primarily regulates immune modulatory processes modulating MS activity. The vitamin D response element was significantly enriched in this system, indicating a direct regulation of this gene interaction network through the vitamin D receptor. With increasing 25(OH)D levels, resulting regulation of this system was associated with a decrease in MS activity. Within the complex network of genes that are regulated by 25(OH)D, well-described targets of IFNB-1b and a regulator of sphingosine-1-phosphate bioavailability were found. The 25(OH)D effects on MS activity were additively enhanced by IFNB-1b.
INTERPRETATION: Here, we provide mechanistic evidence that an unbalanced 25(OH)D gene expression system may affect MS activity. Our findings support a potential benefit of monitoring and managing vitamin D levels (e.g., through supplementation) in early MS patients treated with IFN-beta-1b.

Introduction

There is growing evidence that vitamin D may not only be important in the development of multiple sclerosis (MS), a complex demyelinating disease of the central nervous system, but also in MS disease activity and in its long-term progression.1 Several observational studies of patients with MS have found that increased 25-hydroxyvitamin D (25(OH)D) levels were associated with fewer relapses2–4 and decreased magnetic resonance imaging (MRI) activity.5

Additionally, a number of randomized, placebo-controlled clinical trials of vitamin D supplementation among patients with MS have reported lower numbers of gadolinium-enhancing lesions (GELs)6 and relapses7 or no effect of 25(OH)D on any MS endpoint,8,9 although it must be noted that many clinical trials did not have sufficient power to detect an effect of 25(OH)D on MS disease activity. In a longitudinal study among over 450 participants in the Betaferon/Betaseron in Newly Emerging For Initial Treatment (BENEFIT) clinical trial,10–12 low 25(OH)D levels (<50 nmol/L) predicted conversion from a clinically isolated syndrome (CIS) to clinically definite MS (CDMS) and were associated with a higher rate of new active lesions and a higher relapse rate.13 However, the biological mechanisms underlying these observations are not understood. While 25(OH)D is known to be an immunomodulator,14 how it may act on a molecular level in patients with MS to reduce disease activity is not known.

Thus, we conducted a global gene expression analysis in which expression profiles at various time points among participants in the BENEFIT clinical trial were measured in order to examine the potential relationship between genes or gene sets expressed in association with 25(OH)D and those associated with MS activity. Further analysis of the identified gene sets to characterize their enrichment with respect to known 25(OH)D gene sets and enrichment for vitamin D receptor (VDR) binding sites was performed to help elucidate potential mechanisms for the mode of action of 25(OH)D in the context of MS.

Subjects and Methods

Ethics statement

BENEFIT was a multicenter study, in which patients from 18 countries were included in 98 centers. The study was conducted in agreement with Good Clinical Practice (GCP) principles according to the International Conference on Harmonisation of Technical Requirements for Registration of Pharmaceuticals for Human Use. The Institutional Review Boards of all participating centers approved the study protocol and all patients gave written informed consent before trial entry. The secondary analyses presented in this manuscript were conducted anonymously to describe the impact of vitamin D in the context of interferon beta-1b (IFNB-1b) treatment.

Study population and study design

The BENEFIT study was originally designed to evaluate the impact of early versus delayed treatment with IFNB-1b (Betaseron®; Bayer HealthCare Pharmaceuticals; Berlin, Germany) in patients with a first event suggestive of MS (i.e., a diagnosis of a CIS). Between February 2002 and June 2003, 468 patients from 18 European countries, Israel, and Canada were enrolled at 98 centers. In the initial placebo-controlled phase, patients with a minimum of two clinically silent lesions on MRI, and thus a high probability to convert to CDMS,15 were randomized in a 5:3 ratio to receive either IFNB-1b 250 μg (n = 292) or placebo (n = 176) subcutaneously every other day for 2 years or until diagnosis of CDMS. All patients were then eligible to enter a prospectively planned follow-up phase with open-label IFNB-1b for up to a maximum of 5 years after randomization. Details of the study design and results have been published elsewhere.10–12

Measurement of 25(OH)D

25(OH)D was measured in serum samples collected at baseline (immediately before study treatment in the placebo-controlled phase of BENEFIT and up to 60 days after the onset of the CIS) and at 6, 12, and 24 months thereafter. Under ambient conditions, samples were shipped within 3 days of being drawn and then maintained at −20°C at the central laboratory until further analysis. First thaw process after the initial freeze was completed for this study. Only samples with a minimum of 2 mL of serum were included resulting in 465 patients with at least one, 417 with two or more, 396 with three or more, and 303 with all four 25(OH)D measurements (total of 1570 samples). Serum 25(OH)D was measured using an enzyme-linked immunosorbent assay (ELISA), which is sensitive down to 5.0 nmol/L (Immunodiagnostic Systems Inc., Fountain Hills, AZ). Quality control samples were included among the study samples and the average intraassay coefficient of variation (CV) derived from these samples was 4.4%. The average interassay CV was 11.7%.

Measurement of MS activity

The number of GELs on T1-weighted images from brain MRI scans performed at screening and months 3, 6, 9, 12, 18, 24 and then yearly in the follow-up study was used as a measure of MS activity. There were 1494 simultaneous 25(OH)D and MRI measurements at the CIS (baseline) and months 6, 12, and 24 (Fig. 1). MRI scans were performed as described previously.10,16

Figure 1

Flowchart of available data and size of data partitions used in the different analyses.

Gene expression

At baseline, month 2 or 3, month 12, and month 24 and yearly in the follow-up study, blood for analysis of mRNA expression was collected using the PAXgene system (PreAnalytiX GmbH, Hombrechtikon, Switzerland). Samples were also taken at CDMS if diagnosis occurred before month 24. Gene expression profiles of whole blood were measured using Affymetrix HGU 133 plus 2.0 (affymetrix, Santa Clara, CA, USA) arrays according to the manufacturer's protocol. Raw data were RMA robust multi-array analysis background corrected,17 quantile normalized,18 and transformed with a logarithm to base 2. Probe sets were mapped and summarized to their respective Ensembl gene ID using custom cdf “HGU133Plus2_Hs_ENSG” version 15.1.0 (downloaded from http://brainarray.mbni.med.umich.edu/Brainarray/Database/CustomCDF/15.1.0/ensg.asp, accessed 9 October 2012).19 Global gene expression profiles of ∼19,000 genes were available from 295 individuals; details on available data are shown in Figure 1 and Tables S1, S2.

Statistics

To model the number of GELs as a function of 25(OH)D levels (both as continuous and dichotomized [<50 nmol/L vs. ≥50 nmol/L] variables), IFNB-1b, and other covariates, generalized linear mixed models were applied using proc glimmix in SAS version 9.2 (Cary, NC, USA), which account for the within-person correlation of the repeated measurements. A compound symmetry structure was used to model these correlations. Estimators were exponentiated to give risk ratios (RR). To determine whether 25(OH)D or IFNB-1b was associated with gene expression, generalized linear models were used to model gene expression as a function of 25(OH)D levels, or IFNB-1b, adjusting for age, sex, treatment at the time of 25(OH)D measurement, and the interaction of 25(OH)D level and IFNB-1b. All time points were analyzed together in order to reflect and capture the time-dynamic variation of mRNA and maximize the power of the analysis. The genes that were identified as having expression significantly related to levels of 25(OH)D made up the gene set we named “MS_vitD_genes0.05.” Generalized linear models were also used to model GEL count as a function of gene expression to determine the gene set associated with GEL count. To correct for multiple testing, the Benjamini and Hochberg20 or Bonferroni21 methods were used. For both, a threshold of P ≤ 0.05 was considered statistically significant. R version 2.15.0 was used for all the analyses related to gene expression.

Gene set analysis

Gene Set Enrichment Analysis (GSEA) was done with the JAVA implementation of the GSEA preranked tool, as described previously.22 Briefly, GSEA is a weighted modification of the Kolmogorov–Smirnov test, testing whether the distribution of test scores within a specific gene set significantly differs from the null distribution. Previously published gene sets for comparison were identified from the Molecular Signatures Database version 3.1,22 by searching for genes annotated with the gene ontology (GO) terms GO:0071305 – cellular response to vitamin D; GO:0070563 – negative regulation of VDR signaling pathway; GO:0070564 – positive regulation of VDR signaling pathway; GO:0070562 – regulation of VDR signaling pathway; or GO:0070561 – VDR signaling pathway and from two recently published studies on VDR-mediated gene expression.23,24 A description of each 25(OH)D gene set identified and their respective references are in Table 1. Gene-specific beta estimators quantifying the association of a particular gene with GEL count or the association of 25(OH)D level or IFNB-1b with expression of a particular gene, estimated from the general linear models described above, were used as input vectors in the GSEA tool. For functional gene classification, the Database for Annotation, Visualization, and Integrated Discovery (DAVID)25 was used as follows: the Affymetrix HGU Plus 2.0 platform was selected as background and default DAVID gene sets were tested with the exception of “Functional_Categories” and “Protein_Domains.”

Table 1

Association of 25(OH)D target gene sets with either 25(OH)D, IFNB-1b treatment or Gd-enhancing lesions in BENEFIT, as determined by GSEA

		Association with 25(OH)D		Association with IFNB-1b		Association with GELs
Gene set name	Description	NES	p/FDR-q value1	NES	p/FDR-q value1	NES	p/FDR-q value1
KnightUp [2010]	Genes that are bound and upregulated by the 25(OH)D receptor2	2.09	0.000/0.000	2.46	0.000/0.000	−1.96	0.000/0.004
MS_vitD_genes0.05 Munger et al.	BENEFIT gene expression profiles3	2.08	0.000/0.000	2.21	0.000/0.000	−2.83	0.000/0.000
STAMBOLSKY_RESPONSE_TO_VITAMIN_D3_UP [2010]	Genes induced in SKBR3 cells by 25-hydroxyvitamin D34	1.58	0.006/0.038	1.47	0.039/0.102	1.15	0.236/0.417
CarlbergUp [2011]	Genes that are bound and upregulated by the 25(OH)D receptor5	1.54	0.004/0.040	−1.06	0.357/0.710	0.64	0.974/0.961
BIOCARTA_VDR_PATHWAY [2005]	Control of Gene Expression by Vitamin D Receptor6	1.44	0.088/0.087	−0.93	0.557/0.804	−0.84	0.654/0.983
25(OH)D_GO_associated_genes [2009]	Genes annotated for GO terms related to Vitamin D	1.36	0.121/0.118	1.42	0.102/0.104	0.92	0.581/0.742
V$DR3_Q4 [2005]	Genes with promoter regions containing the motif RRTGNMCYTNNTGAMCCNYNT7	1.01	0.435/0.720	−0.77	0.947/0.844	1.16	0.207/0.527
V$VDR_Q3 [2005]	Genes with promoter regions containing the motif GGGKNARNRRGGWSA8	0.99	0.500/0.672	−1.07	0.267/1.000	1.44	0.002/0.205
V$VDR_Q6 [2005]	Genes with promoter regions containing the motif CNSNNTGAACCN9	0.85	0.828/0.897	0.71	0.978/0.897	1.25	0.067/0.469
STAMBOLSKY_RESPONSE_TO_VITAMIN_D3_DN [2010]	Genes repressed in SKBR3 cells (breast cancer) by 25-hydroxyvitamin D310	0.77	0.814/0.953	1.03	0.418/0.578	1.15	0.301/0.345

25(OH)D, 25-hydroxyvitamin D; IFNB-1b, interferon beta-1b; BENEFIT, Betaferon/Betaseron in Newly Emerging For Initial Treatment; GSEA, Gene Set Enrichment Analysis; FDR, false discovery rate; GELs, gadolinium-enhancing lesions; NES, normalized enrichment score.

FDR < 0.1 designated as statistically significant in GSEA analyses.36

Set derived from Ramagopalan et al.: 25(OH)D receptor ChIP-Seq and GEP in lymphoblastoid cell line after stimulation with calcitriol for 36 h → identification of genes that are bound and upregulated by the 25(OH)D receptor.

Set derived from BENEFIT gene expression profiles and 25(OH)D data from the model GeneExpression∼25(OH)D*IFNß+covariates. To enter the set, the particular gene had to have a nominal P-value of ≤0.05.

Genes induced in SKBR3 cells by 25-hydroxyvitamin D3. From Broad Institute MSigDB v3.1.

Set derived from Heikkinen et al. 25(OH)D receptor ChIP-Seq after 40-min stimulation with calcitriol and gene expression profiling after 4-h stimulation with calcitriol in THP-1 monocytic leukemia cells followed by the identification of genes that are bound and upregulated by the 25(OH)D receptor.

Control of Gene Expression by Vitamin D Receptor. From Broad Institute MSigDB v3.1, derived from Biocarta.

Genes with promoter regions [−2 kb, 2 kb] around transcription start site containing the motif RRTGNMCYTNNTGAMCCNYNT which matches annotation for VDR: vitamin D (1,25- dihydroxyvitamin D3) receptor. From Broad Institute MSigDB v3.1.

Genes with promoter regions [−2 kb, 2 kb] around transcription start site containing the motif GGGKNARNRRGGWSA which matches annotation for VDR: vitamin D (1,25- dihydroxyvitamin D3) receptor. From Broad Institute MSigDB v3.1.

Genes with promoter regions [−2 kb, 2 kb] around transcription start site containing the motif CNSNNTGAACCN which matches annotation for VDR: vitamin D (1,25- dihydroxyvitamin D3) receptor. From Broad Institute MSigDB v3.1.

Genes repressed in SKBR3 cells (breast cancer) by 25-hydroxyvitamin D3. From Broad Institute MSigDB v3.1.

Promoter analysis

Using the Multiple Em for Motif Elicitation (MEME) Suite of motif-based sequence analysis tools,26 two vitamin D-responsive elements (VDRE) positional weighted matrices (PWMs) were inferred from VDR-ChIP sequencing data recently published by Heikkinen et al.23 and Ramagopalan et al.24 (see Fig. S1). Subsequently, we used the Find Individual Motif Occurrences (FIMO) tool with the inferred VDRE PWMs, the VDRE PWM available from the JASPAR database (ID MA0074.1),27 and the consensus Interferon Stimulated Response Element (ISRE) PWM published by Levy et al.,28 to screen for either of these motifs with a P-value detection threshold of 10−4 in the upstream sequences database provided by FIMO. Of note, all three VDRE PWMs are highly similar. However, we screened for all three of them, since this will somewhat account for natural and context-dependent variability of the PWM which is not reflected using point estimates of base frequencies. Using the hypergeometric distribution, we then computed the probability of finding a particular count of VDRE and/or ISRE binding sites in the upstream sequences of the derived 25(OH)D, IFNB-1b and GEL-associated genes given the sites genome-wide background distribution.

Results

Basic characteristics of patients with at least one 25(OH)D measurement are shown in Table 2. 25(OH)D levels by visit are shown in Table 3. The comparison of 25(OH)D levels in patients on IFNB-1b versus placebo at months 6, 12, and 24 did not reveal significant differences between treatment arms. As to be expected, the number of GELs was higher at baseline as compared to the postbaseline visits, with more GELs in the placebo group compared with the IFNB-1b-treated group (Table S3). This difference gradually decreased from months 6 to 24, likely due, in part, to the design of the BENEFIT study which allowed for switching of patients with CDMS in the placebo group to IFNB-1b treatment.

Table 2

Basic characteristics of the cohort of patients from BENEFIT with at least one 25(OH)D measurement

Characteristic	Value
N	465
Age, years ± SD	31.3 ± 7.5
Female sex, n (%)	329 (71)
Number of patients with at least 1 GEL at the CIS, n (%)	196 (42)
Number of patients with at least 1 GEL during the first 24 months, n (%)	262 (56)
IFNB-1b treatment at day 1, n (%)	0 (0)
IFNB-1b treatment at month 24, n (%)	303 (78)

BENEFIT, Betaferon/Betaseron in Newly Emerging For Initial Treatment; 25(OH)D, 25-hydroxyvitamin D; SD, standard deviation; GEL, gadolinium-enhancing lesion; CIS, clinically isolated syndrome; IFNB-1b, interferon beta-1b.

Table 3

Risk for GELs by 25(OH)D level

	N	Mean	Median (IQR)
25(OH)D levels
All	15721	49.53 (22.97)	45.34 (33.43–62.04)
Baseline	457	45.87 (17.93)	43.70 (31.55–56.46)
Month 6	359	48.83 (20.14)	46.65 (34.50–60.64)
Month 12	349	43.38 (20.26)	38.25 (29.17–55.08)
Month 24	386	60.30 (28.71)	54.04 (39.02–76.20)
GEL count2
All	1493	0.89 (2.48)	0.00 (0.00–1.00)
Baseline	453	1.48 (3.38)	0.00 (0.00–1.00)
Month 6	323	0.73 (2.17)	0.00 (0.00–0.00)
Month 12	346	0.66 (1.88)	0.00 (0.00–1.00)
Month 24	371	0.51 (1.68)	0.00 (0.00–0.00)

GEL, gadolinium-enhancing lesion; 25(OH)D, 25-hydroxyvitamin D; IQR, interquartile range; RR, risk ratios.

Note that some 25(OH)D measurements were taken from samples from visits other than baseline, Months 6, 12, or 24. These were excluded from further analysis.

Note that only T1 scans are stated for which a 25(OH)D measurement was available.

Models were adjusted for age at time of measurement, gender, treatment, and the interaction of treatment and 25(OH)D level.

Number of GELs shows inverse association with 25(OH)D levels

An increase of 25(OH)D levels by 50 nmol/L was associated with a statistically significantly lower GEL count (Table 3). Adjusting for gender, age, treatment, and treatment-25(OH)D interaction, the results were unchanged. Similar results were seen when using dichotomized 25(OH)D levels (i.e., groups with 25(OH)D < 50 nmol/L and ≥50 nmol/L, respectively) in the respective model. Significant associations between GEL counts and 25(OH)D levels were also observed when visits were analyzed separately, with the lowest RR for the baseline visit, which had the highest number of GELs and occurred prior to start of IFNB-1b treatment (Table 3). GEL count as a function of 25(OH)D level is shown in Figure S2.

25(OH)D and IFNB-1b systemically regulate 25(OH)D target gene sets

Among BENEFIT participants, the expression of 63 genes was found to be regulated by 25(OH)D based on a nominal P-value threshold of 0.05 (Table S4). However, none of these genes remained statistically significant after false discovery rate (FDR) correction for multiple testing with a P-value threshold of 0.05. These 63 genes made up the 25(OH)D BENEFIT gene set (gene set “MS_vitD_genes0.05” in Table 1). Of note, one of these genes is sphingosine-1-phosphate phosphatase 1 (SGPP1), an enzyme that catalyzes the degradation of sphingosine-1-phosphate (S1P), the substrate for S1P receptors which are a prominent therapeutic target in MS.29,30

The pairwise overlap of the BENEFIT gene set with other 25(OH)D gene sets (see Subjects and Methods and Table 1) was negligibly small (in most cases less than 10%), which supports recent findings that 25(OH)D may regulate a large number of genes and diverse processes in a highly context-dependent manner.24,31–34 By GSEA testing of these 25(OH)D target gene sets, it was determined that the experimentally derived 25(OH)D gene sets (Ramagopalan et al.24 and Stambolsky et al.34) were highly enriched (FDR < 0.05). As all normalized enrichment scores (NES measures the magnitude and direction of a gene set association with the response) were in the positive range, this is indicative of a broad and systemic gene expression upregulation mediated by 25(OH)D (Table 1).

Considering that the rather weak impact of 25(OH)D on single genes (none remained statistically significant after FDR correction) may have been influenced by IFNB-1b-induced gene expression in the BENEFIT cohort, the magnitude of association between all genes and IFNB-1b treatment was examined. This analysis detected 770 genes that were significantly associated with IFNB-1b treatment (Table S5). Sixty-two (98%) of the 63 25(OH)D-associated genes were also significantly associated with IFNB-1b treatment. GSEA analysis of the systemic association of IFNB-1b-induced gene expression with gene sets known to be regulated by 25(OH)D showed that the experimentally derived 25(OH)D target genes sets (MS_vitD_genes0.05, KnightUp,24 and Stambolsky_response_to_vitamin_D3_up34) were highly associated with IFNB-1b gene induction (FDR ≤ 0.1), (Table 1) supporting the hypothesis that a subset of 25(OH)D target genes are coregulated by IFNB-1b in an additive manner.

Genes and gene sets regulated by 25(OH)D are associated with reduction of GELs

Given the association between high 25(OH)D levels with low GEL count, genes regulated by 25(OH)D were examined with respect to their impact on regulating GELs. This could provide a systemic rationale for building a mechanistic model of how 25(OH)D influences MS activity. Modeling GELs as function of gene expression, we found that more than 5000 genes were associated with GEL count (see Table S6 for the top 347 significant genes remaining after Bonferroni adjustment of P-values), and there was overlap with more than 50% of the 25(OH)D-associated genes, as well as 65% of IFNB-1b-associated genes (Fig. 2A). The vast majority of 25(OH)D-associated genes were regulated in such a way by 25(OH)D that would, in turn, reduce GEL count (Fig. 2B–E).

Figure 2

IFNB-1b and 25(OH)D-regulated genes are additionally associated with regulating the count of Gd-enhancing lesions. (A) Venn diagram depicting the overlap between genes that are associated with Gd-enhancing lesions, 25(OH)D or IFNB-1b. (B) Genes whose induction is associated with a reduction of Gd lesion count and which are upregulated through 25(OH)D and IFNB-1b. (C) Genes whose induction is associated with a reduction of Gd lesion count and which are downregulated through 25(OH)D and IFNB-1b. (D) Genes whose induction is associated with an increase of Gd lesion count and which are upregulated through 25(OH)D and IFNB-1b. (E) Genes whose induction is associated with an increase in Gd lesion count and which are downregulated through 25(OH)D and IFNB-1b. IFNB-1b, interferon beta-1b; 25(OH)D, 25-hydroxyvitamin D.

GSEA was performed and revealed that the BENEFIT 25(OH)D gene set (MS_vitD_genes0.05) as well as the KnightUp gene set24 were highly associated with GEL count in such a manner that the upregulation of these gene sets would result in a reduction of GEL count (Table 1).

VDRE and ISRE motifs are enriched in genes associated with 25(OH)D and GELs

Modeling GEL count-related gene expression as a function of 25(OH)D does not provide causal evidence for a relationship and the overlap between 25(OH)D- and IFNB-1b-associated genes may be purely correlative. If there is a causal relationship of these compounds to the genes they putatively regulate, appropriate regulatory elements (namely, VDR response element VDRE and ISRE sites) should be enriched in the promoter regions of the respective genes. Hence, a genome-wide screen with different VDRE and ISRE PWMs using the FIMO tool provided by the MEME suit was conducted.

With respect to the VDRE, these analyses revealed that 8679 of 43267 (Ensembl GRCh37.p10, February 2009 assembly, April 2013 Genebuild) coding and non-coding genes have either of the VDR binding sites present in the sequences upstream of their transcription start site (TSS). This goes well along with the fact that 25(OH)D-mediated VDR activation has been shown to result in a plethora of target genes being regulated.35 Based on this analysis, it was determined that 36 of the 63 genes associated with 25(OH)D in BENEFIT contained at least one VDRE and hence there was a significant enrichment of this regulatory element in the MS_vitD_genes0.05 gene set (P = 0.0012). This indicates that a substantial part of the 25(OH)D-associated genes may directly be regulated by this secosteroid and implies indirect 25(OH)D-mediated regulation of the remaining 27 genes or the presence of regulatory elements in regions that were not included in the analysis (e.g., within exons or introns, the 5′ or 3′ UTR).31 Importantly, 25(OH)D-associated genes in the MS_vitD_genes0.05 gene set also showed significant enrichment for the ISRE motif (P = 0.0009) and additionally the combination of both sites was overrepresented in 25(OH)D-associated genes (P = 0.002). This was in line with the finding that there was a substantial overlap of genes associated with 25(OH)D and genes that were associated with IFNB-1b treatment. For these analyses, FDR values were within the 0.1 cutoff designated as statistically significant in GSEA documentation.36

Next, enrichment for VDRE and ISRE motifs in genes that were significantly associated with GELs was determined. FIMO-based analyses of the upstream sequences of the top genes associated with MS activity revealed that 155 of 347 screened genes contained at least one VDRE motif (P = 0.014) in the sequence upstream of their TSS, 204 contained an ISRE motif (P = 3.8 × 10−17), and 101 contained both a VDRE and ISRE motif (P = 2.5 × 10−6), suggesting a significant enrichment of these sites by themselves and in combination of MS activity-associated genes.

25(OH)D-regulated and GEL-associated genes are functionally enriched for regulation of immune system processes

To complement the single-gene analysis, the biological processes that were regulated by genes that were associated with 25(OH)D in the context of MS activity were inferred. The 63 genes associated with 25(OH)D were analyzed for enrichment of predefined biological terms using DAVID.25 The most prominent enrichment cluster represented gene sets connected to immune responses to pathogens (Table 4). The terms “defense response to bacterium” and “defense response” were detected with nominal P-values of 3.9 × 10−5 and 8.2 × 10−5, respectively. As expected, the top enrichment cluster for IFNB-1b-associated genes also included the term “defense response” (P = 1.1 × 10−10) among other immune system related terms.

Table 4

DAVID enrichment analysis

Term	P-value	Adjusted1 P-value
25(OH)D-associated genes
Defense response to bacterium	<0.001	0.024
Defense response	<0.001	0.026
Response to bacterium	<0.001	0.130
Serine-type endopeptidase activity	<0.001	0.022
Serine-type peptidase activity	<0.001	0.022
IFNB-1b-associated genes
Defense response	<0.001	<0.001
Response to wounding	<0.001	<0.001
Inflammatory response	<0.001	<0.001
Regulation of apoptosis	<0.001	0.004
Regulation of programmed cell death	<0.001	0.005
GEL-associated genes
Defense response	<0.001	0.001
Inflammatory response	<0.001	0.095
MHC class I protein complex	<0.001	0.100
Antigen processing and presentation	<0.001	0.140
Regulation of I-kappaB kinase/NF-kappaB cascade	<0.001	0.140

DAVID, Database for Annotation, Visualization and Integrated Discovery; 25(OH)D, 25-hydroxyvitamin D; IFNB-1b, interferon beta-1b; GEL, gadolinium-enhancing lesions; MHC, major histocompatibility complex.

FDR corrected using the Benjamini and Hochberg (BH) adjustment.20

For GEL count-associated genes, the three most prominent enrichment clusters were all attributable to processes and pathways of the immune system (inflammatory responses, antigen presentation, and NF-kB signaling). The term “defense response” which was also associated with 25(OH)D-associated genes was detected with a P-value of 2 × 10−6. These findings indicate that 25(OH)D regulates processes that are important for the regulation of MS activity.

The 25(OH)D gene interaction network is a component of the MS activity gene interaction network

In order to gain deeper insight into regulatory interactions and associations between the genes found to be associated with GELs and/or 25(OH)D, these two gene sets were analyzed using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING).37 In case of the 25(OH)D-associated genes (Fig. 3A), there is one major connected component centered around lactotransferrin (LTF), an iron-binding protein that plays a major role in immune processes.38 Various genes in this component connected the 25(OH)D system to the IFNB-1b system (e.g., LTF and ELANE). There was a second connected component centered around “interferon-induced protein with tetratricopeptide repeats 1” (IFIT1) which was placed in the central hub-part of the network of genes associated with GELs. This MS network was characterized by a large and densely connected hub component that was particularly enriched for genes connected to interferon response (Fig. 3B).

Figure 3

25(OH)D and MS activity gene interaction networks. The Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) was used to visualize knowledge-based connectivity of genes being associated with 25(OH)D levels (A) or Gd-enhancing lesion count (B). Thickness of lines is proportional to the probability of a true interaction of any pair of genes. 25(OH)D, 25-hydroxyvitamin D; MS, multiple sclerosis.

Discussion

In this large, cross-sectional study of participants in the BENEFIT trial, GEL count was inversely associated with 25(OH)D levels, with every 50 nmol/L increase in 25(OH)D being associated with a 55% decreased GEL rate, although we cannot rule out confounding by other factors that vary with season such as infections. Of note, seasonal variation of GELs and inverse associations with 25(OH)D levels have been previously reported.5,6,39–41

An important gap in the current knowledge is how 25(OH)D may act on a molecular level to influence GEL development. While previous studies have examined the expression of specific 25(OH)D-related genes in patients with MS42,43 and numerous studies on gene expression in MS have been conducted,44–47 this is the first study to scrutinize the relationship between 25(OH)D-driven gene expression and MS activity (as assessed by the number of GELs) over a period of up to 2 years. The results show that even in a cohort such as BENEFIT, with rather low 25(OH)D levels, 25(OH)D exerts systemic gene expression regulation on bona fide 25(OH)D target gene sets derived from lymphocyte and monocyte backgrounds.23 Although on the single gene level, 25(OH)D effects are relatively minor, the results show that this secosteroid is an important systemic factor in that its levels regulate many genes in a subtle fashion which in their sum are important modulators of gene regulatory networks crucial in influencing MS activity.

The analysis of 25(OH)D gene set association with 25(OH)D levels showed strong gene set upregulation, this being suggestive of a context-dependent systemic role of 25(OH)D in reducing GEL count by regulating whole gene sets to jointly influence and tune a variety of biological processes. Of note, this effect appears to be in addition to that of the IFNB-1b treatment-induced gene expression changes, which are known to reduce GEL count10,48 and whose effects have been well described in other studies.44,49,50

Interestingly, some of the 25(OH)D target gene sets are highly significantly associated with reducing GEL count, an effect that is independent of IFNB-1b treatment. This effect is also found when looking at the single genes that are associated with 25(OH)D levels, supporting an additive effect of both 25(OH)D and IFNB-1b in that they reduce GEL count more effectively in combination. Additionally, overrepresentation of VDRE and ISRE motifs in these genes and generally in genes associated with GEL count provides evidence of a mechanistic link. In line with this, the genes associated with 25(OH)D were predominantly associated with immune regulatory processes (in agreement with a recent publication by Satoh and Tabunoki),51 which are also prominently associated with MS and IFNB-1b treatment.44,45,47 One potential mechanism leading to a suppression of MS activity by 25(OH)D could be the S1P receptor-dependent migration of lymphocytes from secondary lymphoid tissues because higher 25(OH)D levels were associated with the expression of SGPP1, a gene that encodes for an enzyme that catalyzes the degradation of S1P. Depriving the S1P receptor of its' normal substrate has been shown to be an effective therapeutic strategy in MS.29,30

25(OH)D involvement in gene expression is known to be context dependent and there are numerous levels at which the mode of action and the magnitude of its effects are fine tuned. Not only does the VDR interact with retinoid X receptor and other cotranscription factors,34 but its expression and the levels of 25(OH)D metabolizing enzymes such as CYP24A1 are also regulated by tightly tuned feedback loops that may be modified by different allelic variants of the molecules involved.52–54 This is also reflected by the quantity of potential 25(OH)D target genes identified in different studies which show only minimal overlap.51 Accordingly, effects of 25(OH)D can be very contrary, reaching from being a proapoptotic agent to being cytoprotective.34 Whereas this study provides evidence that the latter is the case for MS, this ambiguous role of 25(OH)D warrants mechanistic studies of VDR target genes tuning in patients with MS.

As a consequence, these may be used as therapeutic targets that are tunable by 25(OH)D supplementation and that could be a booster of IFNB-1b treatment efficacy.

It is a particular strength of this study that rigorously controlled MRI findings from a clinical trial which was set up to obtain regulatory approval for the use of IFNB-1b in patients with a CIS were used. These data are the first to link disease activity with 25(OH)D-related molecular mechanisms in patients with MS. Because MRI examinations, 25(OH)D measurements, and gene expression profiling were frequently performed simultaneously within the first 2 years of our study, actual inflammatory activity of the disease could be related to 25(OH)D levels and gene expression profiles at the same time point in a cross-sectional manner.

Limitations of this study shall also be addressed. First, expression profiling was performed in whole blood, as the selection of peripheral blood monocytes (or even more specific lymphocyte populations) was not feasible in a multicenter study. This certainly added noise to the findings and it might well be that stronger and more distinct gene expression profiles could have been obtained if more specific cell populations were analyzed. Second, this study focused on gene expression analyses in the cohort treated with IFNB-1b right after CIS; only a limited number of samples obtained under placebo treatment could be included. Therefore, the impact of IFNB-1b on gene expression profiles was largely based on the difference between baseline and postbaseline samples in the early IFNB-1b cohort and, as the baseline sample was obtained right after the CIS, the impact of the first event could be present in those analyses. Nevertheless, the effect of IFNB-1b treatment on the observed gene expression profiles was well in line with published evidence.44,45,47 Third, vitamin D level was only measured at fixed time points (baseline, months 6, 12, and 24), and not at the time of relapses. Relapses at the fixed time points were only observed occasionally. Therefore, the power to compare vitamin-D-related gene expression during relapses and relapse-free periods was too low to perform meaningful analyses. Fourth, although the present analysis suggests a direct effect of 25(OH)D on genes associated with GELs, the possibility of non-VDR-mediated effects that would be associated with increased 25(OH)D and 25(OH)D-associated gene expression, such as ultraviolet light exposure, cannot be excluded and could have independently affected GEL development via non-25(OH)D-related pathways. Finally, almost all patients in this study were whites of European ancestry and the majority of participants were eventually treated with IFNB-1b (although uniform treatment is an important advantage), thus limiting generalizations to individuals of other races or ethnicities or to patients treated with different drugs.

In summary, by analyzing the gene expression, this study identified the first molecular evidence that immunological effects of 25(OH)D impact disease activity in patients with very early MS. These findings should be deepened by the investigation of more specific lymphocyte subsets and extended to populations in the later stages of the disease. Furthermore, biological validation of these results is warranted.