Phenome-wide Mendelian-randomization study of genetically determined vitamin D on multiple health outcomes using the UK Biobank study

Abstract Background Vitamin D deficiency is highly prevalent across the globe. Existing studies suggest that a low vitamin D level is associated with more than 130 outcomes. Exploring the causal role of vitamin D in health outcomes could support or question vitamin D supplementation. Methods We carried out a systematic literature review of previous Mendelian-randomization studies on vitamin D. We then implemented a Mendelian Randomization–Phenome Wide Association Study (MR-PheWAS) analysis on data from 339 256 individuals of White British origin from UK Biobank. We first ran a PheWAS analysis to test the associations between a 25(OH)D polygenic risk score and 920 disease outcomes, and then nine phenotypes (i.e. systolic blood pressure, diastolic blood pressure, risk of hypertension, T2D, ischaemic heart disease, body mass index, depression, non-vertebral fracture and all-cause mortality) that met the pre-defined inclusion criteria for further analysis were examined by multiple MR analytical approaches to explore causality. Results The PheWAS analysis did not identify any health outcome associated with the 25(OH)D polygenic risk score. Although a selection of nine outcomes were reported in previous Mendelian-randomization studies or umbrella reviews to be associated with vitamin D, our MR analysis, with substantial study power (>80% power to detect an association with an odds ratio >1.2 for per standard deviation increase of log-transformed 25[OH]D), was unable to support an interpretation of causal association. Conclusions We investigated the putative causal effects of vitamin D on multiple health outcomes in a White population. We did not support a causal effect on any of the disease outcomes tested. However, we cannot exclude small causal effects or effects on outcomes that we did not have enough power to explore due to the small number of cases.


Introduction
Vitamin D status is an important public-health issue due to the high prevalence of vitamin D insufficiency and deficiency worldwide. 1 Furthermore, it has been reported to be associated with many non-skeletal outcomes (e.g. cardiovascular disease, cognitive impairment and cancer). 2 In our recent umbrella review of meta-analyses of randomized clinical trials (RCTs) and of observational studies, we have found that serum 25-hydroxyvitamin D (25(OH)D) or supplemental vitamin D has been linked to more than 130 unique health outcomes. 3 However, the majority of the studies yielded conflicting results and no association was convincing. 3 With large cohorts linked to electronic medical records (EMRs), the Phenome Wide Association Study (PheWAS) design has been proposed as a high-throughput approach to comprehensively evaluate associations between genetic variants and a wide range of phenotypes (usually generated by EMR). The PheWAS method has been proven to be useful in the replication of hundreds of known genotype-phenotype associations as well the identification of new associations. 4 Since phenotypes defined by EMRs are largely correlated, Bonferroni correction for a conventional PheWAS using general linear models is over-conservative. Thus, a novel Bayesian analysis framework, termed TreeWAS, has been developed. 5 TreeWAS is shown to increase statistical power by up to 20% and can detect new hits missed by a conventional PheWAS. 5 In traditional epidemiological analysis, a causal effect of 25(OH)D on disease outcomes cannot be unambiguously demonstrated due to the known limitations of observational research including unmeasured confounding factors and reverse causality. By using single-nucleotide polymorphisms (SNPs) associated with 25(OH)D in an instrumental variable (IV) analysis, also known as Mendelian-randomization (MR) analysis, we can largely overcome these limitations.
In the current study, we first conducted a systematic literature review of previous MR studies on 25(OH)D. Then, we analysed 339 256 individuals from the UK Biobank study by implementing a conventional PheWAS and the Bayesian TreeWAS for a large number of health outcomes based on linked EMRs. Finally, for a selection of outcomes based on the results from PheWAS/TreeWAS, power and evidence from previous studies, we performed further MR analysis.

Systematic literature review for MR studies on 25(OH)D
The main steps of the study are presented in Figure 1

Study population
UK Biobank is a very large population-based prospective study established to allow detailed investigations of genetic and non-genetic determinants of the diseases in middleand old-aged adults. 6 More than 500 000 individuals were recruited between 2006 and 2010, all of whom gave written consent and underwent baseline measurements, including questionnaires, interview, anthropometric and clinical measurements. Participants donated blood samples for genotyping and biomarker analysis. In addition, UK Biobank participants were linked to their EMR data, including hospital inpatient, cancer-registry and deathregistry data.
In this study, we used the UK Biobank genetic data of 488 378 participants. Genetic quality control was done centrally by UK Biobank. 7

Statistical analysis
We implemented all statistical analyses using R 3.3.2. We used the R package developed by Carroll et al. for the PheWAS analysis. 8   Creation of 25(OH)D genetic-risk score We created a genetic-risk score for 25(OH)D. In the selection of variants, we used the results from the largest genome-wide association study (GWAS) on 25(OH)D in a White population with a total of 79 336 individuals of European ancestry (the SUNLIGHT GWAS). 9 The SUNLIGHT GWAS identified six independent loci associated with serum 25(OH)D concentration that explained 2.84% of the trait variance: rs3755967 (GC), rs12785878 (NADSYN1/DHCR7), rs10741657 (CYP2R1), rs17216707 (CYP24A1), rs10745742 (AMDHD1) and rs8018720 (SEC23A). We created a genetic score by adding the number of effect alleles carried in each of the six SNPs and weighted based on their effect estimates from the SUNLIGHT GWAS.

PheWAS/TreeWAS analysis
The association between the genetic-risk score and common confounders was first tested. In the PheWAS analysis, we only included disease groups with more than 200 cases, as suggested by a simulation of power estimates for PheWAS analysis. 10 We then applied logistic-regression adjusting for gender, age, body mass index (BMI), the UK Biobank assessment centre attended, east and north co-ordinates of home address, and the first five ancestral principal components. A total of 920 disease phenotypes were tested and a P-value of <5.44 Â 10 -5 was regarded as statistically significant based on Bonferroni correction (0.05/920 ¼ 5.44 Â 10 -5 ).
We also applied the TreeWAS Bayesian analysis. The case-control status of participants was defined by their ICD10 codes (from hospital episode, cancer-registry and death-registry data). Due to the complexity of converting ICD9 codes into ICD10 codes, records of ICD9 were discarded. The same weighted genetic-risk score was employed to test the association between the score and all ICD10 codes presented in UK Biobank data. In our TreeWAS analysis, nodes with posterior probability (PP) of a non-zero effect >0.75 were considered as significant. 5

MR analysis
For those phenotypes for which there was enough statistical power for a MR study (>80%) and (i) were statistically significant in the aforementioned PheWAS/TreeWAS analysis or (ii) were classified as probable or suggestive in our previous umbrella review 3 or (iii) were found to be statistically significant or with conflicting evidence in previous MR studies, we implemented MR analysis to further control for bias and assess the causality of the observed associations. To increase the statistical power, we used UK Biobank self-reported medical conditions to include cases that were not captured by EMR data. With power calculation, we estimated that we had 80% power to detect a causal odds ratio (OR) of 1.2 for outcomes with more than 9977 cases assuming that the genetic instrument explains 2.84% of the variance of 25(OH)D levels and that the case:control ratio was 1:5 or larger at alpha ¼ 0.05 level. 11 For these outcomes, we then ran MR analyses using multiple methods: (i) two-stage MR, (ii) inverse variance weighted (IVW) MR and (iii) Egger's regression MR (Supplementary Methods, available as Supplementary data at IJE online). In addition, we constructed the identical genetic-risk score in 2821 control individuals of the Study of Colorectal Cancer in Scotland (SOCCS) (see Supplementary Methods, available as Supplementary data at IJE online, for a description of the SOCCS study) to estimate the variance of 25(OH)D levels explained by this score and the corresponding F-statistics, given that we did not have individual 25(OH)D-level data from the UK Biobank.

Systematic literature review of 25(OH)D MR studies
After applying our inclusion and exclusion criteria, 63 MR studies were included in our systematic literature review (Supplementary Figure 1, available as Supplementary data at IJE online). The causal effect of vitamin D has been examined across a range of disease outcomes and a causal role of vitamin D is not supported for the majority of them. Disease outcomes that were ever reported to be causally associated with vitamin D levels include type 2 diabetes (T2D), total adiponectin, diastolic blood pressure (DBP), risk of hypertension, multiple sclerosis, Alzheimer's disease, all-cause mortality, cancer mortality, mortality excluding cancer and cardiovascular events, ovarian cancer, HDL-cholesterol, triglycerides, high-density lipoprotein, delirium and cognitive functions. However, for some of these outcomes, the evidence across different MR studies is not consistent.

Descriptive analysis of the included UK Biobank participants
We included 339 256 British White unrelated individuals from the UK Biobank cohort, 53.68% of whom were female. All six SNPs satisfied the Hardy-Weinberg equilibrium test (Table 1 and Supplementary Table 3, available as Supplementary data at IJE online). The association between the score of six SNPs and common confounding factors is presented in Table 2. Except for the UK Biobank assessment centre (P ¼ 1.30 Â 10 -17 ), all other confounding factors were not associated with the score.
We then tested the associations between the SNPs and the geographic regions (England and Wales vs Scotland) and found that two SNPs (rs12785878 and rs10745742) were unevenly distributed across the UK. Since we have adjusted for the assessment centre and latitude as covariates, we did not expect the association between score and assessment centre to bias our results.

PheWAS and TreeWAS
In the conventional PheWAS, we tested associations between the score and 920 outcomes (>200 cases). No associations survived the Bonferroni multiple-testing correction. There were only two phenotypes with a P-value <0.001 that were reported as suggestive associations, including delirium (517 cases, P ¼ 1.83 Â 10 -4 ) and nephrotic syndrome (374 cases, P ¼ 9.75 Â 10 -4 ) (Figure 2  and Supplementary Table 4, available as Supplementary data at IJE online). The P-value for the association between the score and vitamin D deficiency was 0.00116 (291 cases), which was the third smallest P-value among all tested associations. We additionally performed a sensitivity analysis without adjustment for BMI (Supplementary Table 5, available as Supplementary data at IJE online) and the results showed no difference in the significance of PheWAS associations when compared to that with adjustment for BMI (Supplementary Table 4, available as Supplementary data at IJE online).
In the TreeWAS, the largest PP for all tested phenotypes was 0.26, which was for the calculus of the ureter (ICD10 code, N20.1). Since a PP >0.75 was needed to take findings from the TreeWAS forward, there were no putative associations found from the TreeWAS analysis.

MR analysis
Based on our eligibility criteria, the following outcomes were selected to be further examined using MR: systolic blood pressure (SBP), DBP, risk of hypertension, risk of T2D, risk of ischaemic heart disease (IHD), BMI, risk of depression, risk of non-vertebral fracture and all-cause mortality. Incorporating self-reported data increased the numbers of cases for some of the outcomes (Table 3). We first applied linear regression between the genetic-risk score and measured plasma 25(OH)D levels in the controls from the SOCCS study (N ¼ 2821). The R 2 value was 1.61% and the F-statistic was 45.96, indicating that the genetic-risk score is a strong IV for our MR analysis. 12 We then applied the second stage of the MR analysis and we did not observe any causal effects for all the tested outcomes. Results from the three different MR methods were consistent (Table 4). Additionally, for the MR Egger's regression, the P-value of the intercept term for all outcomes was >0.05, indicating no evidence of unbalanced pleiotropy (Supplementary Table 6, available as Supplementary data at IJE online) among the variants we used.

Main findings
In order to investigate the safety and rationale of population-wide measurements to raise vitamin D levels (such as fortification of staple food with vitamin D), we conducted a high-throughput PheWAS/TreeWAS study on more than 920 outcomes in a UK Biobank sample of 339 256 British White individuals, followed by MR analyses for selected outcomes. The large sample size ensured the study was well powered to detect moderate to large causal effects (OR >1.2) for outcomes with more than 9997 cases. The PheWAS/TreeWAS analysis used a weighted genetic score of six SNPs as a proxy of the genetically determined 25(OH)D level and examined its association across a wide range of disease outcomes. In this initial scan, we did not observe any significant associations with adjustment for a number of covariates. Additionally, given that there is much debate in the literature suggesting not to adjust for heritable covariates in genome-wide association studies, 13 we performed a sensitivity analysis without BMI as a covariate. There were some differences in the actual effect sizes, but the direction and statistical significance of  the PheWAS associations were consistent between the two models.
We then integrated EMR data with self-reported medical conditions for the final MR analysis and explored the causal effect of vitamin D on SBP, DBP, risk of hypertension, T2D, IHD, BMI, risk of depression, non-vertebral fracture and all-cause mortality. This merged dataset increased the number of cases and re-assigned spurious controls, thus increasing the study power. For outcomes tested in MR, we had statistical power of nearly 100% for detecting a true effect of or larger than 1.2, except for all-cause mortality (85% power). The estimated causal-effect sizes for all outcomes were close to null with narrow confidence intervals, which suggested that there is no moderate to large causal effect of 25(OH)D on the nine tested outcomes.
Among these nine outcomes, associations between vitamin D and blood pressure, hypertension, T2D, IHD and BMI were suggestive from previous observational studies and RCTs. 3 Conflicting evidence exists for blood pressure and T2D from previous MR studies. [14][15][16][17][18] Previous MR studies for IHD and BMI did not support causal associations. [19][20][21][22] Although associations between vitamin D and the risk of depression and non-vertebral fracture were suggestive from observational studies and RCTs, our MR analysis found no evidence of causality for their associations. These finding are further supported by the recently published MR studies, 23,24 in which the effects of vitamin D on major depression (59 851 cases) 23 and fracture (185 057 cases) 24 are investigated, respectively, but none of them supports any causality. In agreement with MR findings, a recent systematic literature review and meta-analysis investigating the effect of calcium, vitamin D or a combination of calcium and vitamin D supplements on the incidence of fractures (33 RCTs included, involving 51 145 participants) did not find any association between vitamin D or calcium plus vitamin D supplements and the incidence of non-vertebral fractures. 25 There was no evidence on the association between vitamin D and all-cause mortality from previous observational studies and RCTs, although a previous MR study reported a significant effect [10 349 deaths, 95 766 total participants, OR ¼ 1.30, 95% confidence interval (CI): 1.05 to 1.61]. 26 This study employed four SNPs in two loci (DHCR7 and CYP2R1) as their IV, which explained only 1.0% of the variance for the 25(OH)D level. However, in our study with 85% power for detecting an effect of 1.2 and with a comparable case size (N ¼ 9830 deaths), we did not observe a causal effect (OR ¼ 1.030, 95% CI: 0.869 to 1.222, P ¼ 0.671). The association between vitamin D and all-cause mortality needs to be studied by a larger MR or meta-analysis of MR studies.

Strengths and limitations
Our study has several strengths. This is the first study to investigate the causal effect of vitamin D in a large sample of 339 256 individuals across the whole spectrum of disease outcomes. Taking advantage of the PheWAS design, we tested the association between the 25(OH)D genetic-risk score and a wide spectrum of phenotypes. Then, we applied multiple MR methods, including two-stage MR, IVW MR and Egger's MR, to explore the causal effect and test the robustness of our findings across multiple methods.
The study also has some limitations. In our analysis, we only included White individuals residing at a high latitude, which may hinder the generalizability of our findings to other populations. The weights we employed in the score creation were from a meta-analysis of GWAS, 9 which covered individuals residing in multiple countries. Therefore, the distribution of vitamin D levels in the UK population may be different from that in the wider White population and thus the actual coefficients for variants might also differ. Another implication of the high latitude is that, in the UK Biobank population, an overall low 25(OH)D level (independent of genetic variation but due to a lack of adequate sunlight exposure) may increase the risk of vitamin D-related disease. Therefore, small changes in 25(OH)D levels due to genetic variation may not affect the disease risks further. Furthermore, the variance of 25(OH)D explained by the six SNPs was 2.84% from the SUNLIGHT GWAS and 1.61% from the SOCCS controls. Although the F-statistics indicates a robust instrument, the variance explained is low. Moreover, considering participants in the UK Biobank are not representative of the UK population and there is evidence of a 'healthy volunteer' selection bias, 27 it is possible that the healthy volunteer selection bias observed may contribute to the null findings in this study. Furthermore, the relationship between the 25(OH)D level and the risk of diseases may be nonlinear. As shown by previous studies, vitamin D supplementation only shows treatment effects among individuals with baseline 25(OH)D levels of no more than 30 nmol/L. When all participants were analysed irrespective of their baseline 25(OH)D levels, there was no treatment effect. 28,29 Thus, the effect of 25(OH)D on health outcomes may differ by baseline serum 25(OH)D level. Considering the potential divergent 25(OH)D levels of the UK population, it is possible that we missed the true association between 25(OH)D levels and diseases among individuals of certain 25(OH)D levels.
Although we incorporated EMR data and self-reported medical conditions in our definition of phenotypes, problems with reporting bias could have occurred in outcome definition. Incorporation of more data, including general-practice data, outpatient data, prescription data and even imaging data, would help improve the validity of case definition in PheWAS studies. We set the minimum case number per phenotype based on a simulation analysis of PheWAS power estimates. 10 We therefore restricted the PheWAS analysis to outcomes with >200 cases. From PheWAS analysis, we noted that the association between the genetic score of the 25(OH)D levels and vitamin D deficiency was not significant. This is probably due to the limited statistical power. In particular, only 291 cases with vitamin D deficiency were identified in this cohort of 339 256 individuals, which is much fewer than would be expected if 25(OH)D levels were systematically tested. Additionally, for the MR analyses, we explored only nine outcomes with more than 80% power. There were some outcomes that had been shown to be associated with vitamin D in previous MR studies, including total adiponectin, 30 multiple sclerosis, 31 Alzheimer's disease, 32 cancer mortality, 26 mortality excluding cancer and cardiovascular events, 26 ovarian cancer, 33 HDL-cholesterol, 16 triglycerides 16 and cognitive functions, 34 but, due to limited statistical power or data availability, we did not include them in our final MR analysis.

Conclusion
Our study suggested that there was no evidence of a large to moderate (OR >1.2) causal effect of vitamin D on a number of health outcomes, particularly for SBP, DBP, the risk of hypertension, T2D, IHD, BMI, depression, nonvertebral fracture and all-cause mortality. Further, larger studies, probably involving the joint analysis of data from several large biobanks, may be needed to investigate smaller causal effects that nevertheless could be important for public health due to the high prevalence of low vitamin D levels in many populations.