Mendelian randomization study updates the effect of 25-hydroxyvitamin D levels on the risk of multiple sclerosis

Observational studies and previous Mendelian randomization (MR) studies have shown that genetically low 25-hydroxyvitamin D (25OHD) levels are associated with a high susceptibility to multiple sclerosis (MS). The present MR study aims to update the causal estimates for the effects of 25OHD levels on MS risk. To date, the largest genome-wide association study (GWAS) for serum 25OHD (n = 401,460) and MS (14,498 MS cases and 24,091 controls) was used to assess the effect of serum 25OHD levels on MS. All participants were of European ancestry. The MR-egger_intercept test and Cochran’s Q statistic were used to determine the pleiotropy and the heterogeneity, respectively. MR-egger, weighted median, inverse variance weighted (multiplicative random effects), simple mode, and weighted mode methods were used to evaluate the causal association of serum 25OHD levels with MS. Finally, the effect of a single 25OHD SNP (single nucleotide polymorphism) on MS was used to test the SNP bias. One hundred and fifteen newly identified serum 25OHD genetic variants were extracted from a large-scale serum 25OHD GWAS dataset. The 20 most effective and independent 25OHD genetic instrumental variables were extracted from the MS GWAS summary statistics. Pleiotropy analysis suggested no significant pleiotropic variant among the 20 selected 25OHD genetic instrument variants in MS GWAS datasets. As serum levels of 25OHD based on genetic changes increased, the risk of MS decreased using MR-egger (Beta = − 0.940, p = 0.001; OR = 0.391), weighted median (Beta = − 0.835, p = 0.000; OR = 0.434), IVW (Beta = − 0.781, p = 0.000; OR = 0.458), simple mode (Beta = − 1.484, p = 0.016; OR = 0.227), and weighted mode (Beta = − 0.913, p = 0.000; OR = 0.401). Our results were robust, with no obvious bias based on investigating the single 25OHD SNP on MS. Our analysis suggested a causal association between genetically increased serum 25OHD levels and reduced MS in the European population.

pathogenic mechanisms driving MS [5]. It is known that MS is mediated by the interaction between genetic and environmental factors and well-powered genome-wide association studies (GWAS) have been explored to investigate the genetic risk factors in MS [6].
The active form of vitamin D, 1,25-Dihydroxyvitamin D3 (1,25(OH)2D3), suppresses autoimmune diseases such as MS by reducing the production of proinflammatory cytokines such as IFN (interferon)-γ, IL (interleukin)-2, and IL-17, as well as enhancing the secretion of anti-inflammatory cytokines such as IL-4 and IL-10 [7][8][9][10][11]. A recent retrospective case-control study showed that the effect of vitamin D3 on MS susceptibility was not mediated by regulatory T cells (Tregs) [12]. Increased evidences suggest that decreased vitamin D levels is associated with an abnormal immune response in MS [13,14]. Low vitamin D levels have been shown to correlate with disease activity in various autoimmune disorders and neurological diseases including MS [13][14][15][16][17][18][19]. Numerous epidemiological studies have strongly suggested that vitamin D insufficiency contributes to MS risk [20]. One SNP, rs10766197 CYP2R1 which regulates serum 25-OHvitamin D3 (25OHD) levels, is reported to relate with MS risk [21,22]. Thus, hypovitaminosis D was considered as a possible risk factor for MS, with multiple decades of research on the association of vitamin D status and MS [23].
Mendelian randomization (MR) studies have been widely used to assess the causal link between an exposure and an outcome by an analytical method using genetic variants as instrumental variables [13,[24][25][26]. Four previous MR studies used three or four serum 25OHD genetic instrumental variables (IVs) and suggested a link between predicted serum 25OHD and the risk of MS [27][28][29][30]. MR studies have many limitations, including low statistical power. Thus, to improve statistical power, 20 serum 25OHD genetic IVs from the largest GWAS for serum 25OHD (n = 401,460) were used to update the effect of serum 25OHD levels on MS risk.

Ethics approval and consent to participate
Our study was approved by the Ethics Committee of Beijing Institute of Brain Disorders in Capital Medical University. This article contains human participants collected by several studies to report the large-scale GWAS for the serum 25OHD [31] and MS [32]. All participants gave informed consent in all the corresponding original studies, as described in the Methods.

Study design
Three principal assumptions were prerequisites in the MR study [24,25] and are described in Fig. 1. Assumption 1 is that genetic instrument variants are reliably associated with serum 25OHD levels. Assumption 2 is that genetic instrument variants should undoubtedly not be associated with any confounders. Assumption 3 is that genetic instrument variants firmly influence the risk of MS through the 25OHD and not through other pathways.

Genetic instrument variants for serum 25OHD
Sixty-nine independent 25OHD genetic variant loci were identified using imputed genotypes from 401,460 white British UK Biobank participants with available serum 25OHD levels [31]. The TwoSampleMR R package was used to perform a two-sample MR study. One hundred and fifteen independent 25OHD genetic variants with each trait at p < 5 × 10 −8 are demonstrated in Additional file 1: Table S1 and were selected as "suggestive" instruments in our MR studies, as described in recent studies [24,25]. To obtain a more accurate effect size estimate with these variants, the following rules was used to Fig. 1 Design flow chart for the MR study. MR assumptions: assumption 1, 2 and 3. Solid line represents direct causal effects that 25OHD genetic instrument variants are reliably associated with 25OHD levels and influence the risk of MS through the 25OHD in assumption 1. The dotted line represents that 25OHD genetic instrument variants are not associated with any measured or unmeasured confounders and do not influence the risk of MS through other pathways in assumption 2 and 3, respectively age, smoking, and socioeconomic position) were also removed. Finally, 20 effective and independent 25OHDassociated single nucleotide polymorphisms (SNPs) were selected as the potential IVs and are shown in Table 2.

MS GWAS dataset
The MS GWAS dataset is from the International Multiple Sclerosis Genetics Consortium (IMSGC) and consists of 14,498 MS subject and 24,091 healthy control individuals of European ancestry [32]. The TwoSampleMR R package was used to perform a two-sample MR study. Twenty effective and independent 25OHD genetic variants were extracted from MS GWAS dataset (https:// gwas. mrcieu. ac. uk/ datas ets/ ieu-a-1025/) by using functions "extract_ outcome_data" and "harmonise_data". Table 1 provides the sample profiles for the MS GWAS datasets.

Pleiotropy test
The MR-egger_intercept test was used to determine the pleiotropy test [33]. The TwoSampleMR R package was used to perform the MR-egger_intercept test. The MR-egger_intercept test was performed using function "mr_pleiotropy_test". P < 0.05 is the level of statistical significance.

Heterogeneity test
Cochran's Q statistic was used to determine the heterogeneity test [34,35]. The TwoSampleMR R package was used to perform Cochran's Q statistic. The heterogeneity test was performed by using function "mr_heterogeneity". P < 0.05 is the level of statistical significance.

MR analysis
The TwoSampleMR R package was used to perform the MR analysis. The MR analysis was performed by using  function "mr". We selected five MR analysis methods including MR-egger, weighted median, inverse variance weighted (multiplicative random effects), simple mode, weighted mode [33,36,37]. The effect size (beta) and 95% confidence interval (CI) correspond to 1 standard deviation (SD) in 25OHD levels. P < 0.05 is the level of statistical significance.

Single SNP effect analysis
The TwoSampleMR R package was used to analyze the effect of serum single 25OHD SNP on MS. This analysis was performed by using functions "mr_singlesnp" and "mr_leaveoneout". Forest plot were prepared by using functions "mr_forest_plot" and "mr_leaveoneout_plot".

Serum 25OHD genetic instrument variants in MG GWAS dataset
To explore the effect of serum 25OHD genetic instrument variants on MS, we used the MS GWAS dataset (https:// gwas. mrcieu. ac. uk/ datas ets/ ieu-a-1025/), described in Table 1, to extract 20 effective and independent 25OHD genetic variants from the MS GWAS dataset by using functions "extract_outcome_data" and "harmonise_data". Table 2 showed the detailed information about 20 25OHD genetic instrument variants in the MS GWAS dataset.

Pleiotropy analysis
MR-egger_intercept test was performed to determine the pleiotropy test by using function "mr_pleiotropy_test". Table 3 shows detailed information about the pleiotropy analysis of 20 25OHD genetic instrument variants in MG GWAS dataset. The statistical P = 0.409 suggests that there is no significant pleiotropic variant among the selected 20 25OHD genetic instrument variants in the MS GWAS datasets. Thus, we established the assumption that genetic instrument variants should not be associated with any confounders. Therefore, all these 20 25OHD genetic variants could be used as the effective instrumental variables in this MR analysis.

Heterogeneity test
Cochran's Q statistic was performed to test the heterogeneity by using function "mr_heterogeneity". Table 4 shows detailed information about the heterogeneity analysis of 20 25OHD genetic instrument variants in the MS GWAS dataset. Both MR-egger and inverse variance weighted methods demonstrated the statistical P = 0.04. P < 0.05 suggests significant heterogeneity among the 20 25OHD genetic instrument variants in MG GWAS datasets. Therefore, we mainly use the random effect model to estimate the MR effect size.

MR analysis
We found that as 25OHD levels increased due to genetic differences, the risk of MS decreased. We  Table 5). The individual MR estimates demonstrate that as the effect of a single SNP on 25OHD increased, the suppressive effect of a single SNP on MS increased, as determined using five MR analysis methods including MR-egger, weighted median, inverse variance weighted (multiplicative random effects), simple mode, and weighted mode (Fig. 2). Collectively, our data suggest causal association of increased 25OHD levels from greater genetic variation with decreased risk of MS.

Single SNP effect analysis
A single effect of serum 25OHD SNP was performed by using functions "mr_singlesnp" and "mr_leaveoneout". Forest plots were prepared by using functions "mr_for-est_plot" and "mr_leaveoneout_plot". A forest plot of the 20 serum 25OHD SNPs associated with risk of MS is demonstrated in Fig. 3. MR leave-one-out sensitivity analysis for the effect of the 20 serum 25OHD SNPs on MS is shown in Fig. 4. These findings indicate that our results are robust and have no obvious bias based on the effect of the single 25OHD SNP on MS.

Discussion
Epidemiological evidence suggests an association between 25OHD deficiency and risk and disease progression of MS [13,20,[38][39][40][41][42][43][44]. This relationship is further proven by four previous MR studies using three or four serum 25OHD genetic IVs [27][28][29][30]. In the first MR study in 2015, Mokry et al. [28] identified an odds ratio [OR] of 2 by using four serum 25OHD genetic IVs obtained from a serum 25OHD GWAS (N = 33,996) [45], as well as a MS GWAS involving 14,498 MS cases and 24,091 healthy controls from International Multiple Sclerosis Genetic Consortium (IMSGC) [32]. In the second MR study in 2016, Rhead et al. [29] identified an OR = 0.79 using three serum 25OHD genetic IVs. These samples were obtained from a serum 25OHD GWAS (N = 4501) [46] and MS GWAS from members of Kaiser Permanente Medical Care Plan, Northern California Region (KPNC) (1,056 MS cases and 9,015 controls) and the Epidemiological Investigation of Multiple Sclerosis and Genes and Environment in Multiple Sclerosis MS case-control studies (6335 cases and 5762 controls). In the third MR study in 2020, Jacobs et al. [27] identified an OR = 0.57 using six serum 25OHD genetic IVs obtained from a serum 25OHD GWAS (N = 79,366) [47] and IMSGC discovery phase GWAS (14,802 MS,26,703 controls). In the fourth MR study in 2021, Harroud et al. [30] identified an OR = 0.72 using four serum 25OHD genetic IVs obtained from a serum 25OHD GWAS (N = 33,996) and IMSGC discovery phase GWAS (14,802 MS, 26,703 controls). In the present study, we identified an OR = 0.22 ~ 0.45 using 20 serum 25OHD genetic IVs from to date, the largest serum 25OHD GWAS (n = 401,460) [31] and MS GWAS (14,498 MS cases and 24,091 controls) [32]. Thus, our MR analysis further strengthens evidence for a causal link between 25OHD levels and MS risk.
Among 115 newly identified serum 25OHD genetic variants, 20 effectively independent 25OHD genetic instrumental variables were extracted from MS GWAS summary statistics. We did not identify any significant pleiotropic variant among the selected 20 25OHD genetic instrument variants in the MS GWAS datasets by pleiotropy analysis. Of note, five MR analysis methods demonstrated that 25OHD genetic instrument variants firmly influence the risk of MS through 25OHD but not through other pathways. We found that increased 25OHD levels due to genetic changes were significantly associated with reduced MS risk. For every 1 SD, genetically increased 25OHD levels could reduce MS risk 0.940, 0.835, 0.781, 1.484, and 0.913 using MR-egger, weighted median, inverse variance weighted (multiplicative random effects), simple mode Table 5 The causal association of serum 25OHD levels with MS IVW Inverse variance weighted. IVW was based on multiplicative random effects. nsnp the number of single-nucleotide polymorphism. Beta: the regression coefficient based on 25OHD raising effect allele. SE: standard error. p < 0.05 represents the causal association of serum 25OHD levels with MS. OR Odds ratio. OR_lci95 Lower limit of 95% confidence interval for OR. OR_uci95 Upper limit of 95% confidence interval for OR  and weighted mode, respectively. These results suggest a causal association between genetically increased 25OHD levels and MS risk. Both MR-egger and inverse variance weighted methods demonstrated a significant heterogeneity among the 20 25OHD genetic instrument variants in the MS GWAS datasets. Thus, we mainly used the random effect model to estimate the MR effect size. Inverse variance weighted (multiplicative random effects) showed a reduced trend of MS with 1 SD genetically increased 25OHD levels (beta = − 0.781, 95% CI: [− 1.105, − 0.457], p = 2.34E−06). MR-egger, weighted median, simple mode and weighted mode also proved a causal link between 25OHD levels and MS risk. In addition, a single SNP effect analysis and leave-one-out sensitivity analysis for the effect of the 20 serum 25OHD SNPs on MS suggests 25OHD deficiency is becoming an increasing problem worldwide. A recent finding led to a proposed link between 25OHD deficiency and autoimmune diseases [48]. It is reported that 25OHD deficiency is associated with MS risk [39]. Apart from the maintenance of healthy calcium metabolism, Vitamin D and its active product (1,25-dihydroxyvitamin D) have additional roles in immune and central nervous system cell homeostasis [38].
To date, there is evidence to suggest there are some benefits of vitamin D for patient-important outcomes among people with MS [40]. However, vitamin D did not affect recurrence of relapse and worsening of disability in patients with MS [40]. However, in clinical trials, while doses of up to 40,000 IU/day were tested and appeared safe as add-on therapy for short-term periods, chronic and high-dose therapy can lead to life-threatening complications related to vitamin D toxicity, including renal failure, cardiac arrythmia and status epilepticus [43]. There were several strengths to this MR study. First, this MR design was based on three principal assumptions and we used different methods to prove three principal assumptions. Second, we used the large-scale serum 25OHD GWAS dataset (n = 401,460) and MS GWAS summary statistics (14,498 MS cases and 24,091 controls) [31,32]. With a large-scale GWAS, it is easy to find a causal association between 25OHD levels and MS risk. Third, the two GWAS from European ancestry reduce the influence of population stratification. Fourth, using independent statistical methods, we proved 20 independent genetic variants to be the effective instrumental variables. Fifth, we used five MR analysis methods including MRegger, weighted median, inverse variance weighted (multiplicative random effects), simple mode, and weighted mode. Critically, all five MR methods proved the causal association between 25OHD levels and MS risk. Lastly, our results were robust, with no obvious bias based on the different methods.
Of course, our MR study had some limitations. First, our MR analysis is only based on European ancestry. Therefore, our conclusions need to be investigated in other ancestries. Second, the GWAS dataset from IMSGC is based on clinically diagnosed MS [32]. Different diagnostic criteria may produce different results. Third, among 115 newly identified serum 25OHD genetic variants [31], only 20 25OHD genetic instrumental variables were effective and independent. Fourth, although we proved that increased serum 25OHD levels could reduce MS risk, the mechanisms by which 25OHD reduced MS risk is still unclear. Fifth, the genetically increased serum 25OHD level does not necessarily reflect the active form of vitamin D, 1,25-Dihydroxyvitamin D3 (1,25(OH)2D3), in the serum [49]. Sixth, while our MR studies may provide the strongest evidence for a causal link between vitamin D levels and the risk of MS, our study does not address the role of vitamin D on disease activity in patients with established MS. Therefore, our study must expand to explore the causal association of vitamin D levels with the activity of MS disease using other databases.
Although MR study offers several advantages, it also has many limitations [50]. These limitations include pleiotropy, low statistical power, canalization, population stratification (which is confounding by ancestry), etc. [50]. To detect and correct for pleiotropy, probably the most challenging to address, a robust method MR-egger intercept was used in our MR study. To improve the statistical power, we used the, to date, largest large-scale GWAS for serum 25OHD [31] and MS [32]. To eliminate population stratification, both 25OHD GWAS and MS GWAS were chosen from European descent. Although we try to minimize the shortcomings, maybe there are still some limitations. Thus, our results need be reproduced in randomized controlled trials (RCTs).

Conclusions
In summary, our MR analysis strengthens evidence for the causal association between genetically increased serum 25OHD levels and reduced MS in European people.