Genetic Prediction of Serum 25-Hydroxyvitamin D, Calcium, and Parathyroid Hormone Levels in Relation to Development of Type 2 Diabetes: A Mendelian Randomization Study

Vitamin D and parathyroid hormone (PTH) are the two major regulators of serum calcium (S-Ca) (1). Without vitamin D, only 10–15% of dietary calcium is absorbed, and vitamin D reduces renal calcium losses. PTH mobilizes calcium from the skeleton when S-Ca is reduced, and PTH also improves renal conservation of calcium. In addition to these S-Ca–regulating effects, vitamin D and PTH, and S-Ca itself, have been suggested to influence the development of several different diseases.

The association between serum 25-hydroxyvitamin D (S-25OHD) and the risk of type 2 diabetes (T2DM) is unclear. A large-scale meta-analysis including 21 observational prospective cohort studies (76,220 participants and 4,996 incident cases) with average follow-up periods from 1.3 to 22 years revealed a significant inverse association between S-25OHD levels and risk of T2DM across different populations (2). However, the T2DM-protective effect of vitamin D was not observed in randomized controlled trials (3,4). In contrast to S-25OHD, observational studies have reported an increased risk of T2DM associated with elevated S-Ca (5) and S-PTH levels (6), but whether these associations are causal is uncertain due to the potential of residual confounding and reverse causality.

Using genetic variants as instrumental variables for an exposure, Mendelian randomization (MR) studies can strengthen the inference on the causal nature of exposure-outcome relationships (7). There are several outstanding strengths of the MR design compared with standard observational studies (7). According to Mendel’s Law of Inheritance, genetic variants are randomly assorted at conception and therefore less likely to be associated with potential confounding factors. It also diminishes reverse causality, since genotype cannot be affected by the development of disease. In addition, genetically predicted exposure reflects lifelong exposure to the potential risk factor. There are three basic assumptions of an MR study (7) (Supplementary Fig. 1). First, the genetic variants selected as instrumental variables should be robustly associated with the exposure. Second, the instrumental variables should not be associated with any confounders. Third, the genetic variants that predict the exposure should affect the risk of the outcome merely through the exposure, not via other alternative pathways.

Findings of previous MR studies of S-25OHD levels instrumented by two or four single nucleotide polymorphisms (SNPs) in relation to T2DM are inconclusive, with an inverse association observed in some studies (8,9) but no association in others (10,11). Recently, two novel S-25OHD–related SNPs were detected in a large genome-wide association study (GWAS) (12). To clarify the associations of S-25OHD and S-Ca levels with T2DM, we conducted an updated two-sample MR study, using a larger sample size for T2DM and more instrumental variables compared with previous MR studies of S-25OHD levels in relation to T2DM (8–11). In order to explore the potential mechanisms behind the associations, the associations of S-25OHD and S-Ca levels with four glycemic traits (fasting glucose, fasting insulin, and HOMA of β-cell function [HOMA-B] and insulin resistance [HOMA-IR]) were assessed in secondary analyses. In addition, we explored the association between S-PTH levels and T2DM, which has, to the best of our knowledge, not been examined previously using the MR design.

This two-sample MR study is based on summary-level data from the published GWAS on S-25OHD, S-Ca, and S-PTH levels and T2DM. Detailed information of the GWAS and SNPs used as instrumental variables is presented in Supplementary Tables 1–3. The analyses have been approved by the Swedish Ethical Review Authority.

Six SNPs for S-25OHD, seven SNPs for S-Ca, and five SNPs for S-PTH were selected at the genome-wide significance level (P < 5 × 10⁻⁸) from the summary-level data of the hitherto largest GWAS on these traits, which included 79,366, 39,400, and 29,155 individuals of European ancestry, respectively (12–14). These SNPs explained ∼5.3%, 0.9%, and 4.5% of the variation in S-25OHD, S-Ca, and S-PTH levels, respectively (12–14). For each trait, the included SNPs were located in different gene regions and distributed independently (i.e., not in linkage disequilibrium). For all three traits, adjustments had been made for age, sex, and study-specific covariates if needed (e.g., principal components of ancestry, study site, and season) (12–14). S-25OHD levels had further been adjusted for BMI and assay batch (12). The vitamin D assay used differed among the included studies (12). For S-25OHD levels, an additional low-frequency SNP (rs10741657 in CYP2R1) of large effect on S-25OHD levels was included from a GWAS of 42,274 individuals of European ancestry (15), resulting in seven instrumental variables for S-25OHD. No linkage disequilibrium was detected between this low-frequency SNP with another S-25OHD–related SNP in CYP2R1 (pairwise R² = 0.03) (15). Three of the SNPs associated with S-25OHD levels are located in or near genes (CYP2R1 and DHCR7) involved in vitamin D synthesis, whereas the other four SNPs are located in or near genes involved in the transport (GC and SEC23A) or catabolism (CYP24A1 and amidohydrolase domain–containing 1 [AMDHD1]) of vitamin D. For S-Ca, the SNP in the GCKR gene (Supplementary Table 2) was removed due to multiple pleiotropic associations with potential confounders (16), leaving six SNPs as instrumental variables for S-Ca levels.

Summary-level data for the association of the SNPs influencing S-25OHD, S-Ca, and S-PTH levels with T2DM were available in the DIAbetes Genetics Replication And Meta-analysis (DIAGRAM) consortium, which includes 32 studies with a total of 898,130 individuals (74,124 case subjects with T2DM and 824,006 control subjects) of European ancestry (17) (Supplementary Table 1). Genetic association data for T2DM adjusted for BMI were used in the primary analyses, and data without BMI adjustment were used in sensitivity analyses. Summary-level data for fasting glucose, fasting insulin, HOMA-B, and HOMA-IR were obtained from the Meta-Analyses of Glucose and Insulin-related traits Consortium (MAGIC) (18). Data on fasting glucose were based on a meta-analysis of 21 GWAS including 46,186 individuals of European ancestry without diabetes (18). Summary-level data of fasting insulin, HOMA-B, and HOMA-IR were acquired from meta-analysis of 20 GWAS with, respectively, 38,238, 36,466, and 37,037 individuals of European descent without diabetes (18). One SNP (rs117913124) associated with S-25OHD was not found in the data sets of glycemic traits, and no suitable proxy SNP was available.

The random-effects inverse variance–weighted method was used for the main analyses. Sensitivity analyses applying the fixed-effects inverse variance–weighted, weighted median (19), and MR-Egger regression methods (20) were performed. The PhenoScanner database V2 was searched to investigate potential pleiotropic associations of the instrumental variables (21). The odds ratios (ORs) of T2DM were scaled to 1 SD increments in genetically predicted S-25OHD, S-Ca, and S-PTH levels. For glycemic traits, β estimates were calculated per SD increment of plasma S-25OHD and S-Ca levels. An approximate SD for S-25OHD and S-PTH levels was obtained from the population-based Swedish Mammography Cohort and corresponded to 0.33 ln-nmol/L for S-25OHD and 0.33 ln-pg/mL for S-PTH (22). The SD for calcium corresponded to 0.5 mg/dL (equivalent to 0.125 mmol/L) and was obtained from the GWAS for S-Ca (13).

In order to investigate the effect of S-25OHD synthesis on T2DM, three SNPs in or nearby genes related to S-25OHD synthesis (rs117913124, rs12785878, and rs10741657) were used as instrumental variables (8). Associations with P values <0.017 (where P = 0.05/three exposures) were deemed statistically significant after Bonferroni correction for three exposures. P values between 0.017 and 0.05 were regarded as suggestive evidence of associations. All P values were two sided. The statistical analyses were performed in Stata/SE 15.0.

For a genetically predicted 1 SD change in serum levels, we had 100% statistical power to detect an OR of 0.9 (or 1.1) for T2DM in the analyses of S-25OHD and S-PTH and an OR of 0.8 (or 1.2) in the analyses of S-Ca. However, we only had ∼70% power to detect an OR of 0.9 (or 1.1) for T2DM in the analyses of S-Ca.

There was suggestive evidence of an inverse association of genetically predicted S-25OHD levels with T2DM in an analysis based on all seven SNPs and the random-effects inverse variance–weighted method. For 1 SD increment of S-25OHD levels, the OR of T2DM was 0.94 (95% CI 0.88–0.99; P = 0.029), with moderate heterogeneity between SNPs (I² = 38%; P = 0.138) (Fig. 1). Results from the weighted median method supported an inverse association, but the results were not statistically significant (Fig. 1). The MR-Egger regression analysis revealed no association and provided evidence of directional pleiotropy (P = 0.012) (Fig. 1). In a search of the PhenoScanner database, the SNPs rs3755967 and rs17216707 were found to be associated with immune system and kidney function, respectively (Supplementary Table 4). rs10741657 was associated with height, hip circumference, and appendicular muscle mass (Supplementary Table 4). No other significant pleiotropic effects were found. An analysis based on BMI-unadjusted data showed a nonsignificant inverse association between S-25OHD levels and T2DM (OR 0.96 [95% CI 0.92–1.01]; P = 0.088) (Supplementary Fig. 2).

Restricting the analysis to the three S-25OHD–related SNPs within or near genes in the vitamin D synthesis pathway revealed an OR of T2DM of 0.90 (95% CI 0.83–0.98; P = 0.011) for each SD increase in S-25OHD levels, without heterogeneity among SNPs (I² = 5%; P = 0.351) (Fig. 2). The results were supported in an analysis based on the weighted median method (Fig. 2). The MR-Egger regression analysis showed no significant association, but the precision was low and there was no evidence of pleiotropy (P = 0.153). An analysis based on BMI-unadjusted T2DM data yielded similar results (OR 0.93 [95% CI 0.88–1.00]; P = 0.04).

There was no association between genetically predicted S-Ca levels, based on six SNPs, and T2DM in the random-effects inverse variance–weighted analysis (OR 1.09 [95% CI 0.73–1.65]; P = 0.665) (Fig. 3). Heterogeneity was observed between the estimates of individual SNPs (I² = 65%; P = 0.013) (Fig. 3). Among SNPs, rs1801725 in CASR was positively associated with T2DM, whereas rs7481584 in CARS was inversely associated with T2DM (Fig. 3). The overall null association was found in the weighted median model (Fig. 3). The MR-Egger regression analysis showed evidence of directional pleiotropy (P = 0.002) and yielded an OR of 1.41 (95% CI 1.12–1.77; P = 0.003) per 1 SD increment of S-Ca levels (Fig. 3). The association remained positive after exclusion of rs1801725 in CASR (OR 2.12 [95% CI 0.43–10.56]; P = 0.357) and after further exclusion of rs7481584 in CARS (OR 1.63 [95% CI 0.32–8.31]; P = 0.554), but results were not statistically significant. Results from the PhenoScanner search revealed possible pleiotropic factors related to SNPs of S-Ca levels (Supplementary Table 5). Specifically, the calcium-raising allele of rs1801725 in the CASR gene region was associated with lower diastolic blood pressure and pulse rate. However, diastolic blood pressure has been found to be positively associated with risk of T2DM (23) and therefore cannot explain or mediate the association with rs1801725. Rs1550532 and rs7481584 were associated with bilirubin levels and inflammatory response, respectively (Supplementary Table 5), which are related to T2DM (24,25).

No association between genetically predicted S-PTH levels and T2DM was observed, but there was substantial heterogeneity between the five SNPs (Fig. 4). Among SNPs, rs73186030 (highly correlated with the S-Ca–associated SNP rs1801725; R² = 0.96; P < 0.0001) in the CASR gene was significantly associated with higher odds of T2DM (Fig. 4). SNPs in CYP24A1 and RGS14 were also positively and inversely, respectively, associated with T2DM (P < 0.05) (Fig. 4).

No associations were found between genetically predicted S-25OHD and S-Ca levels and four glycemic traits (Supplementary Table 2 and Supplementary Figs. 3 and 4). Only two of the five PTH-associated SNPs were available in the glycemic trait data sets. Because of this and the lack of association between S-PTH and T2DM, no analyses of the association of S-PTH levels with the glycemic traits were performed.

Genetically predicted S-25OHD levels were inversely associated with T2DM in this MR study based on 898,130 European-descent individuals from 32 studies. The inverse association was robust in the analysis confined to the three SNPs within or close to genes in the vitamin D synthesis pathway. There was evidence of horizontal pleiotropy in the analysis based on all S-25OHD–related SNPs. This study found evidence of a positive association between genetically predicted S-Ca levels and T2DM after adjusting for horizontal pleiotropy. The exploratory analysis of S-PTH levels in relation to T2DM showed no overall association.

Epidemiological studies on the association between S-25OHD and T2DM have provided conflicting results (2,10–13). Some individual observational studies (2,10,11) and several meta-analyses of observational studies (2,5) have reported an inverse association between circulating vitamin D levels and risk of T2DM in different populations. However, randomized controlled trials have found no significant effect of vitamin D supplementation alone (3,26,27) or vitamin D plus calcium supplementation (4) on incident T2DM over a median follow-up time ranging from 1 to 7 years among healthy individuals or populations at high risk of T2DM (e.g., prediabetes) (3,26), osteoporotic fracture (27), and hypovitaminosis D (26). One of the trials found a suggestive protective effect, with a hazard ratio of 0.88 (95% CI 0.75–1.04; P = 0.12), comparing vitamin D supplementation alone with placebo (3). A possible explanation for the lack of significant association is the relatively short follow-up, especially in one of the trials (3), and low power. Moreover, the addition of calcium in one of the trials (4) may have attenuated any potential protective effect of vitamin D supplementation. Other potential limitations are that some trials were conducted in populations at high risk of T2DM (prediabetes) (3,26) or osteoporotic fracture (27) and that nonadherence in placebo or intervention groups may have shifted the risk difference between groups toward the null. To the best of our knowledge, none of the trials tested for effect modification by SNPs influencing S-25OHD levels. The Vitamin D and Omega-3 Trial (VITAL), a large trial involving 25,871 U.S. adults, has been conducted to examine the effects of vitamin D (D3) supplements on T2DM prevention. Findings of T2DM-related outcomes have not yet been published, but recent results from the trial showed that vitamin D supplementation did not significantly reduce the risk of major cardiovascular events or cancer but reduced the risk of cancer-related deaths (28).

Findings of previous MR studies are also inconsistent (8–11). An individual-level MR analysis including five prospective cohorts (n = 28,144 cases) found no association between S-25OHD levels instrumented by four SNPs and incidence of T2DM in European populations (11). A null association was also found in an Asian population (10). However, in line with the results of our study, other MR analyses based on data from three Danish cohorts (n = 5,037 case subjects) and the China Kadoorie Biobank (n = 5,565 case subjects) as well as a meta-analysis of the China Kadoorie Biobank and nine studies of European individuals (n = 58,312 case subjects in total) revealed a significant inverse association of genetically predicted S-25OHD levels instrumented by two SNPs in the vitamin D synthesis pathways with T2DM (8,9).

As for potential mechanisms, the current study found no support that the association between S-25OHD levels and T2DM may be mediated by fasting glucose, fasting insulin, islet β-cell function, or insulin resistance. These results are in agreement with most experimental and longitudinal studies showing no benefit of vitamin D supplementation on improving hyperglycemia, β-cell secretion, or insulin sensitivity in healthy or high-risk populations or in patients with T2DM with or without S-25OHD deficiency (29–31).

Other potential mechanisms whereby high levels of S-25OHD might reduce T2DM risk are through upregulation of certain immune pathways or inhibition of inflammation (32,33). A systemic review of cell studies revealed that elevated S-25OHD levels suppressed certain vital marker expression, such as macrophage chemotactic protein 1, interleukin 6, and interleukin 8, and exerted anti-inflammatory effects (32). Reviews based on clinical studies have also found that S-25OHD levels may change the balance between inflammatory cytokines, thereby inhibiting inflammation and decreasing T2DM risk (33). The activated form of vitamin D (1,25-dihydroxyvitamin D₃) may influence the function of pancreatic cells through facilitating the expression of insulin receptor and enhancing insulin responsiveness for glucose transport (34). S-25OHD can regulate calcium fluxes in the β-cells, thereby regulating insulin secretion (34).

Among the seven SNPs used as instrumental variables for S-25OHD levels in this MR study, those within or near genes related to vitamin D synthesis (DHCR7 and CYP2R1) and catabolism (CYP24A1) as well as a SNP in the AMDHD1 gene showed evidence of an inverse association with T2DM (Fig. 1). DHCR7 encodes the enzyme 7-dehydrocholesterol reductase, which converts 7-dehydrocholesterol to cholesterol, thereby removing the substrate from the synthesis pathway of vitamin D₃, a precursor of 25OHD (35). CYP2R1 encodes 25-hydroxylase, the key enzyme in the conversion of vitamin D₃ to 25OHD. CYP24A1 encodes 24-hydroxylase, which catalyzes reactions involved in degradation of both 25OHD and 1,25-dihydroxyvitamin D, the active form of vitamin D (35). AMDHD1 encodes an enzyme involved in the catabolic pathway of histidine, lysine, phenylalanine, tyrosine, proline, and tryptophan (12). Although most of the S-25OHD–related SNPs that were associated with T2DM are located within or close to genes related to vitamin D metabolism, we cannot rule out that the observed associations between genetically predicted S-25OHD levels and T2DM are driven by pleiotropic effects related to, for example, cholesterol and amino acids.

Studies on the association of S-Ca levels with T2DM are limited and inconsistent (5,6,36). The Atherosclerosis Risk in Communities study included 12,800 U.S. adults and observed a multivariable-adjusted hazard ratio of T2DM of 1.26 (95% CI 1.07–1.48; P = 0.004), comparing the highest and lowest S-Ca level quintile (6). Elevated S-Ca levels were also associated with an increased risk of diabetes in a retrospective cohort study of 6,096 Hong Kong Chinese adults (5). However, in a 10-year cohort study including 8,800 Korean adults, albumin-adjusted S-Ca levels were not associated with T2DM risk after multivariable adjustment (36). In the current study, genetic evidence supports a possible association of S-Ca levels with T2DM, but the association was only observed after adjustment for pleiotropy. In addition, among the six instrumental variables for S-Ca levels, only two showed a positive association with T2DM, of which the SNP in the calcium-sensing receptor (rs1801725 in CASR) was significantly positively associated with T2DM. A positive association between that SNP and T2DM has also been shown in a previous study (6). A GWAS including up to 22,653 participants revealed the importance of rs73186030 (highly correlated with rs1801725) in the CASR gene in regulating S-PTH and S-Ca levels (16). Whether the association between genetic variations in CASR and T2DM is related to S-Ca or S-PTH levels and their downstream effects, or to horizontal pleiotropy bias, remains unclear.

A major strength of our study is the MR study design, which minimized residual confounding and reverse causality by using genetic variants as proxies of lifelong higher S-25OHD, S-Ca, and S-PTH levels. In addition, population stratification was diminished because the analyses were based on data from individuals of European ancestry, which confines the transferability of the present findings to other populations. Thus, studies of a causal nature are warranted to verify the present conclusions among individuals of different ancestries. We used summary-level data from large-scale genetic consortia, thereby assuring high statistical power to detect weak associations. A limitation of this study is that we found evidence of directional pleiotropy in the MR-Egger regression analysis for both S-25OHD (when using all seven associated SNPs) and S-Ca levels. However, the S-25OHD–associated SNPs did not seem to have pleiotropic associations with potential confounders. Although two of the S-Ca–associated SNPs were associated with possible confounders, the association remained positive, though nonsignificant, after excluding those SNPs. A potential shortcoming is that the study population came from the general population. Our findings may therefore not be generalizable to at-risk populations or individuals with vitamin D deficiency. Another limitation is that we were unable to examine whether there are any interactions between genetically predicted S-25OHD, S-Ca, and S-PTH levels and lifestyle/environmental factors on the risk of T2DM.

The current study provides some genetic support that S-25OHD levels are inversely associated with risk of T2DM. However, a consistent association between S-25OHD levels and T2DM was only observed in an analysis based on SNPs in the vitamin D synthesis pathway. High S-Ca levels, primarily caused by genetic variations in the calcium-sensing receptor, may be associated with an increased risk of T2DM.

Acknowledgments. Summary-level data for genetic associations with T2DM were contributed by the DIAGRAM consortium. Summary-level data for genetic associations with fasting glucose, fasting insulin, HOMA-B, and HOMA-IR were contributed by the MAGIC consortium. The authors thank all investigators for sharing these data.

Funding. This work was supported by the Swedish Research Council and the Research Council for Health, Work Life, and Welfare.

Duality of Interest. No potential conflicts of interest relevant to this article were reported.

Author Contributions. S.Y. wrote the manuscript. S.Y., X.J., K.M., and S.C.L. interpreted the data and reviewed the manuscript. S.Y. and S.C.L. designed the study and analyzed the data. S.C.L. is the guarantor of this work and, as such, had full access to all the data in the study and takes responsibility for the integrity of the data and the accuracy of the data analysis.