Coronavirus disease 2019 (COVID-19), which is caused by the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), has become a global pandemic^[1]. Primary prevention, which is aimed at delaying disease onset by modulating modifiable risk factors, has been proposed to address the global challenge posed by COVID-19.

Vitamin D plays a critical role in the immune system, and vitamin D deficiency is a risk factor for several adverse health outcomes, such as cancer, cardiovascular diseases, metabolic disorders, infectious diseases, and autoimmune diseases. Some studies have shown that low 25 hydroxyvitamin D (25OHD) concentrations are associated with a high risk of COVID-19^[2], but other studies have found no association^[3]. Therefore, the association between 25OHD concentrations and the risk of COVID-19 has not been established. The above studies were observational, and traditional epidemiological studies are vulnerable to reverse causality and residual confounding.

A promising approach, known as Mendelian randomization (MR), which uses inherited genetic variants as variables, can provide better evidence for the causal effect of exposure on diseases and largely overcome the traditional limitations related to confounding and reverse causality^[4]. The two-sample MR method is an excellent strategy for evaluating causality using summary statistics from genome-wide association study (GWAS) data. The summary statistics from the GWAS data released by the UK biobank and Severe COVID-19 GWAS Group is based on the European population. The data were obtained from the largest sample, as of now, and they were aimed for using in the exploration of the genetic variants of the outcomes related to COVID-19 susceptibility and severity^[5]. In the present study, we conducted an MR analysis to investigate a causal relationship between 25OHD concentration and COVID-19.

The summary data were obtained from the most recent version of the results of the GWAS analyses for the COVID-19 host genetics initiative involving UKB individuals, which was released on July 1, 2020 (https://www.covid19hg.org/results/). Two phenotypes, including COVID-19 and severe respiratory infection confirmed as COVID-19, were examined. Detailed information on the phenotypes is provided in Supplementary Table S1 (available in www.besjournal.com). The summary data on COVID-19 and severe respiratory infection confirmed as COVID-19 were based on three comparisons: patients with COVID-19 (n = 6,696) and the general population without the phenotype (n = 1,073,072)—the population hereafter was named UKB-COVID-19; patients with COVID-19 (n = 3,523) and those without COVID-19 (n = 36,634)—the population hereafter was named UKB-COVID-19-negative; patients with severe respiratory infection confirmed as COVID-19 (n = 536) and the general population without the phenotype (n = 329,391)—the population hereafter was named severe UKB-COVID-19. In addition, we used summary data for severe respiratory confirmed COVID-19 reported by the Severe COVID-19 GWAS Group (https://ikmb.shinyapps.io/COVID-19_GWAS_Browser/). The study included 1,160 patients with severe respiratory infection confirmed as COVID-19 and 2,205 control participants from the general population without COVID-19 (hereafter severe COVID-19). The definitions of severe respiratory infection confirmed as COVID-19 by the Severe COVID-19 GWAS Group and the COVID-19 host genetics initiative were different (Supplementary Table S1).

One hundred and forty-three single-nucleotide polymorphisms (SNPs) were selected as instrumental variables (IVs) for the 25OHD concentration based on a recent large-scale GWAS involving 417,580 Europeans^[6], which explains the 10.5% variance of the serum 25OHD concentrations. In addition, we retained independent variants (linkage disequilibrium [LD], r² < 0.001) for sensitivity analysis. The SNP with the lowest P-value was selected when we encountered genetic variants with LD. The LD proxies were defined using 1,000 European genomes.

For the analyses, only the genetic variants with available SNP-exposure and SNP-outcome association data were retained. The inverse-variance weighted (IVW) method was used to estimate the overall association between 25OHD concentrations and COVID-19 susceptibility and severity, and the weighted median, penalized weighted median, and MR-Egger regression were used to detect potential violations of the valid instrumental variable assumptions. The MR-Egger analysis was performed to evaluate pleiotropy based on the intercept. The heterogeneity tests were performed using IVW Q-test. Sensitivity analyses were performed to examine the stability of the causal estimate. Firstly, we performed a “leave one out” analysis to further investigate the possibility that the association was driven by a single SNP. We retained independent variants (LD, r² < 0.001) for further sensitivity analysis. The results are presented as odds ratios (OR) with 95% confidence interval (CI) and beta with standard error (SE) of outcomes per genetically predicted increase in each exposure factor.

Regarding the various estimates for different measures, we chose the result of the main MR method according to the following rules:

(1) If there was no directional pleiotropy in the MR estimates (Q statistic: P-value > 0.05, MR-Egger intercept: P-value > 0.05), the results of the IVW method were reported.

(2) If directional pleiotropy was detected (MR-Egger intercept: P-value < 0.05) and a P-value of > 0.05 was observed for the Q-test, the results of the MR-Egger method were reported.

(3) If directional pleiotropy was detected (MR-Egger intercept: P-value < 0.05) and a P-value of < 0.05 was observed for the Q-test, the results of the weighted median method were reported.

All data analyses were performed using the “TwoSampleMR” package for R version 4.0.0 (https://www.r-project.org/).

As shown in Table 1, the MR analysis showed no significant association between the genetically modified 25OHD concentration and COVID-19 in the population of UKB-COVID-19 and UKB-COVID-19 patients (OR = 1.136, 95% CI: 0.988−1.306, P = 0.074; OR = 1.168, 95% CI: 0.956−1.427, P = 0.128). The association between 25OHD concentrations and COVID-19 was robust based on the outcomes of the weighted median and penalized weighted median methods, except for the MR-Egger regression (OR = 1.258, 95% CI: 1.053−1.502, P = 0.013; OR = 1.302, 95% CI: 1.011−1.676, P = 0.044). Pleiotropy bias and heterogeneity were not observed. We report the results of the IVW method based on various estimates for different measures. In addition, the “leave one out” results showed that no individual genetic variants seemed to have any significant effect on the overall results by omitting the included 89 SNPs one at a time (Supplementary Figures S1−S2 available in www.besjournal.com). In addition, the non-causal association between 25OHD concentrations and COVID-19 remained robust based on assessment with the SNP instrument (LD, r² < 0.001) (Table 1, Supplementary Figures S3−S4 available in www.besjournal.com).

The MR analysis showed no significant association between the genetically modified 25OHD concentrations and severe COVID-19 in the population with severe UKB-COVID-19 and severe COVID-19 (OR = 0.889, 95% CI: 0.549−1.439, P = 0.246; OR = 0.894, 95% CI: 0.587−1.363, P = 0.603) (Table 2). This emphasized that the weighted median, penalized weighted median, and MR-Egger regression methods were robust. There was limited evidence of heterogeneity and horizontal pleiotropy based on the Q-test and MR-Egger intercept test. Sensitivity analyses using different instruments yielded similar findings, suggesting the robustness of the observed insignificant association (Table 2, Supplementary Figures S5−S8 available in www.besjournal.com).

The present two-sample MR analysis demonstrated that 25OHD concentrations did not appear to be associated with the risk of COVID-19 susceptibility and severity. Four previous similar MR analyses were performed to explore the association between 25OHD concentrations and the risk of COVID-19^[7-10]. Of these, four used the summary GWAS data for vitamin D from the SUNLIGHT Consortium, three also used the summary GWAS data for vitamin D from UKB^{[7, 9, 10]}, and two used individual data^{[7, 10]}. Compared with the four previous studies, our study used the GWAS summary data for the 25OHD concentrations with the largest explained variance (10.5%) from the UKB and found a non-causal association between 25OHD and these verity of COVID-19 based on the summary data from the Severe Covid-19 GWAS Group. Our findings were consistent with those of three of the above studies^[8-10]; they used individual data ^[10] and performed MR analysis using summary GWAS data for vitamin D from the SUNLIGHT Consortium.

There are some possible explanations for these negative findings. First, only a few of the genetic variants that affect 25OHD concentrations are included, and the findings may have been biased due to weak IVs. Second, these null findings suggest that the association between 25OHD concentrations and COVID-19 susceptibility and severity may be attributed to reverse causation bias and confounder bias. Vitamin D from the environment, through sunlight or diet, is metabolized in the liver to 25OHD, which is used to determine the vitamin D status of a patient^{[6, 11]}. Vitamin D deficiency is common in patients with severe COVID-19^[2], which may be a consequence of quarantine and reduced outdoor behavior or acute infection in COVID-19 patients. The above findings suggest that COVID-19 susceptibility and severity are expected to decrease the prevalence of vitamin D deficiency, which needs to be proven by more bidirectional MR studies. Third, it is also important to note that the MR study considers the lifelong effect of the genetic modification of COVID-19. However, the association between vitamin D concentrations and the risk of COVID-19 may vary with time. The cross-sectional observational nature of all the current MR studies limits the evaluation. Future MR studies incorporating follow-up data should consider the effect of vitamin D concentrations on COVID-19 and how changes in the genetic variant effects with time may impact the interpretability and validity of their results.

This study has some limitations. The MR analyses were based on numerous assumptions. First, the genetic variants were associated with 25OHD concentrations. We selected genetic variants as IVs, based on a recent large-scale GWAS^[6], which showed strong associations with 25OHD concentrations; therefore, the bias associated with weak instruments may be less likely. Second, the genetic variants were not associated with measured and unmeasured confounders that influenced both vitamin D and COVID-19. However, the unmeasured confounders (i.e., actual vitamin D intake and UVB exposure) or alternative causal pathways (i.e., the expression of vitamin D binding protein) may still have affected our results because of the limitation of the method. Third, the horizontal pleiotropy may distort the MR results. In our study, there was limited evidence of heterogeneity and horizontal pleiotropy. Fourth, all associations are linear. However, it is difficult to validate this assumption. The violation of this assumption is not essential when the aim is to test the null hypothesis of no effect of vitamin D concentration on COVID-19, but it can cause problems when the aim is to assess the size of the effect. In addition, the GWAS of severe COVID-19 cases involved a small sample, which may have had an effect on the MR estimate and limited the IVs for COVID-19 for reverse MR analysis. There is an overlap between the UK Biobank and the COVID-19 Host Genetics Institute GWAS samples. However, our study adopted a two-sample MR analysis, which may have resulted in a bias. The two samples were selected from the European population and a mixed population, which may also have resulted in bias. In addition, the limitations of relationship between the changes in vitamin D concentrations over time and COVID-19 risk should be considered. Therefore, future studies with larger samples and individual data in one sample are needed to verify and explore the observed associations.

In conclusion, we performed a two-sample MR analysis to explore the associations between vitamin D and the risk of COVID-19 to address the limitations of confounding and reverse causality in observational studies. Our findings showed that 25OHD concentrations did not appear to be associated with the risk of COVID-19 susceptibility and severity. In the future, bi-directional MR should be performed to further validate these results. Future MR studies incorporating follow-up data should consider the effect of vitamin D concentrations on COVID-19 and how changes in genetic variants with time may impact the interpretability and validity of their results.

Figure S1.  MR leave-one-out sensitivity analysis for ‘25OHD concentration’ on ‘COVID-19’ in the population of UKB-COVID-19 based on 89 SNPs. Leave-one-out analysis: each row represents a MR analysis of 25OHD concentration on COVID-19 using all instruments expect for the SNP listed on the y-axis. The point represents the beta with that SNP removed and the line represents 95% confidence interval. The summary data are reported by the COVID-19 host genetics initiative. COVID-19, coronavirus disease 2019; MR, mendelian randomization; SNP, single-nucleotide polymorphism; 25OHD, 25-hydroxyitamin D; SD, standard deviation; UKB, UK biobank

Figure S2.  MR leave-one-out sensitivity analysis for ‘25OHD concentration’ on ‘COVID-19’ in the population of UKB-COVID-19 negative based on 89 SNPs. Leave-one-out analysis: each row represents a MR analysis of 25OHD concentration on COVID-19 using all instruments expect for the SNP listed on the y-axis. The point represents the beta with that SNP removed and the line represents 95% confidence interval. The summary data are reported by the COVID-19 host genetics initiative. COVID-19, coronavirus disease 2019; MR, mendelian randomization; SNP, single-nucleotide polymorphism; 25OHD, 25-hydroxyitamin D; SD, standard deviation; UKB, UK biobank.

Figure S3.  MR leave-one-out sensitivity analysis for ‘25OHD concentration’ on ‘COVID-19’ in the population of UKB-COVID-19 based on 56 SNPs. Leave-one-out analysis: each row represents a MR analysis of 25OHD concentration on COVID-19 using all instruments expect for the SNP listed on the y-axis. The point represents the beta with that SNP removed and the line represents 95% confidence interval. The summary data are reported by the COVID-19 host genetics initiative. COVID-19, coronavirus disease 2019; MR, mendelian randomization; SNP, single-nucleotide polymorphism; 25OHD, 25-hydroxyitamin D; SD, standard deviation; UKB, UK biobank

Figure S4.  MR leave-one-out sensitivity analysis for ‘25OHD concentration’ on ‘COVID-19’ in the population of UKB-COVID-19 negative based on 57 SNPs. Leave-one-out analysis: each row represents a MR analysis of 25OHD concentration on COVID-19 using all instruments expect for the SNP listed on the y-axis. The point represents the beta with that SNP removed and the line represents 95% confidence interval. The summary data are reported by the COVID-19 host genetics initiative. COVID-19, coronavirus disease 2019; MR, mendelian randomization; SNP, single-nucleotide polymorphism; 25OHD, 25-hydroxyitamin D; SD, standard deviation; UKB, UK biobank

Figure S5.  MR leave-one-out sensitivity analysis for ‘25OHD concentration’ on ‘severe COVID-19’ in the population of severe UKB-COVID-19 based on 82 SNPs. Leave-one-out analysis: each row represents a MR analysis of 25OHD concentration on severe COVID-19 using all instruments expect for the SNP listed on the y-axis. The point represents the beta with that SNP removed and the line represents 95% confidence interval. The summary data are reported by the COVID-19 host genetics initiative. COVID-19, coronavirus disease 2019; MR, mendelian randomization; SNP, single-nucleotide polymorphism; 25OHD, 25-hydroxyitamin D; SD, standard deviation; UKB, UK biobank

Figure S6.  MR leave-one-out sensitivity analysis for ‘25OHD concentration’ on ‘severe COVID-19’ in the population of severe COVID-19 based on 82 SNPs. Leave-one-out analysis: each row represents a MR analysis of 25OHD concentration on severe COVID-19 using all instruments expect for the SNP listed on the y-axis. The point represents the beta with that SNP removed and the line represents 95% confidence interval. The summary data are reported by the severe COVID-19 GWAS Group. COVID-19, coronavirus disease 2019; MR, mendelian randomization; SNP, single-nucleotide polymorphism; 25OHD, 25-hydroxyitamin D; SD, standard deviation; GWAS, genome wide association study

Figure S7.  MR leave-one-out sensitivity analysis for ‘25OHD concentration’ on ‘severe COVID-19’ in the population of severe UKB-COVID-19 based on 55 SNPs. Leave-one-out analysis: each row represents a MR analysis of 25OHD concentration on severe COVID-19 using all instruments expect for the SNP listed on the y-axis. The point represents the beta with that SNP removed and the line represents 95% confidence interval. The summary data are reported by the COVID-19 host genetics initiative. COVID-19, coronavirus disease 2019; MR, mendelian randomization; SNP, single-nucleotide polymorphism; 25OHD, 25-hydroxyitamin D; SD, standard deviation; UKB, UK biobank

Figure S8.  MR leave-one-out sensitivity analysis for ‘25OHD concentration’ on ‘severe COVID-19’ in the population of severe COVID-19 based on 53 SNPs. Leave-one-out analysis: each row represents a MR analysis of 25OHD concentration on severe COVID-19 using all instruments expect for the SNP listed on the y-axis. The point represents the beta with that SNP removed and the line represents 95% confidence interval. The summary data are reported by the severe COVID-19 GWAS Group. COVID-19, coronavirus disease 2019; MR, mendelian randomization; SNP, single-nucleotide polymorphism; 25OHD, 25-hydroxyitamin D; SD, standard deviation; GWAS, genome wide association study