A Bidirectional Mendelian Randomization Study to evaluate the causal role of reduced blood vitamin D levels with type 2 diabetes risk in South Asians and Europeans

Context Multiple observational studies have reported an inverse relationship between 25-hydroxyvitamin D concentrations (25(OH)D) and type 2 diabetes (T2D). However, the results of short- and long-term interventional trials concerning the relationship between 25(OH)D and T2D risk have been inconsistent. Objectives and methods To evaluate the causal role of reduced blood 25(OH)D in T2D, here we have performed a bidirectional Mendelian randomization study using 59,890 individuals (5,862 T2D cases and 54,028 controls) from European and Asian Indian ancestries. We used six known SNPs, including three T2D SNPs and three vitamin D pathway SNPs, as a genetic instrument to evaluate the causality and direction of the association between T2D and circulating 25(OH)D concentration. Results Results of the combined meta-analysis of eight participating studies showed that a composite score of three T2D SNPs would significantly increase T2D risk by an odds ratio (OR) of 1.24, p = 1.82 × 10–32; Z score 11.86, which, however, had no significant association with 25(OH)D status (Beta -0.02nmol/L ± SE 0.01nmol/L; p = 0.83; Z score -0.21). Likewise, the genetically instrumented composite score of 25(OH)D lowering alleles significantly decreased 25(OH)D concentrations (-2.1nmol/L ± SE 0.1nmol/L, p = 7.92 × 10–78; Z score -18.68) but was not associated with increased risk for T2D (OR 1.00, p = 0.12; Z score 1.54). However, using 25(OH)D synthesis SNP (DHCR7; rs12785878) as an individual genetic instrument, a per allele reduction of 25(OH)D concentration (-4.2nmol/L ± SE 0.3nmol/L) was predicted to increase T2D risk by 5%, p = 0.004; Z score 2.84. This effect, however, was not seen in other 25(OH)D SNPs (GC rs2282679, CYP2R1 rs12794714) when used as an individual instrument. Conclusion Our new data on this bidirectional Mendelian randomization study suggests that genetically instrumented T2D risk does not cause changes in 25(OH)D levels. However, genetically regulated 25(OH)D deficiency due to vitamin D synthesis gene (DHCR7) may influence the risk of T2D. Supplementary Information The online version contains supplementary material available at 10.1186/s12937-021-00725-1.


Introduction
Type 2 diabetes mellitus (T2D) has become a global health epidemic of twenty-first century. The International Diabetes Federation (IDF) data showed that the number of people affected with T2D would rise from 382 million in 2011 to 592 million by 2035 [1]. T2D is a complicated disease impacted by the complex interplay of genetic, epigenetic, and environmental factors [2][3][4]. In addition, several lifestyle factors such as sedentary lifestyle, westernized diet, smoking, and genetic predisposition enhance T2D risk in some ethnic groups [5,6].
Vitamin D deficiency increases in conjunction with T2D, Type 1 diabetes mellitus (T1D), obesity, and cardiovascular disease [7]. The inverse relationship between circulating 25(OH)D concentrations and T2D is confirmed in several observational and prospective studies from the US, Australia, Europe, and Asia [8][9][10][11]. However, a causal link has not been established between low blood 25(OH) D levels and T2D risk. Recently, a large clinical trial of vitamin D3 supplementation on 2423 participants for 2.5 years concluded that a dose of 4000 IU of 25(OH)D3 did not protect people from developing T2D [12]. Findings from observational studies may be subject to residual confounding because it is difficult to measure individual variation specifically related to sun exposure and cultural variations in such studies. It may also be challenging to rule out the reverse causality accounting for any association. On the other hand, human genetic information is used by the methodology of Mendelian randomization, which takes account of genetic instruments to provide an unconfounded estimate of the association. The Mendelian randomization strategy is based on the principle that the individual genotypes are randomly assigned and can be used as a genetic instrument with the assumption that their involvement affects the outcome only by modifying the biomarkers (i.e., circulating 25(OH)D) and can also be used to test the direction of causation [13,14].
The blood level of 25(OH)D is the best and commonly used method to determine vitamin D status. According to the US Endocrinology Society, vitamin D deficiency is defined as a 25(OH)D level of ≤ 20 ng/mL; insufficiency as 21 to 29 ng/mL; and sufficiency as > 30 ng/mL or more [15]. Low 25(OH)D status (deficiency or insufficiency) is highly prevalent among elderly, postmenopausal women, and certain ethnic population groups and is influenced by both genetic and environmental factors [16][17][18][19]. Mainly, populations due to specific ethnic backgrounds (South Asian, East Asian, and African) show a high prevalence of blood 25(OH)D deficiency [10,20,21]. Genetically regulated 25(OH)D concentration in a longitudinal follow-up study of 95,766 Europeans of Danish ancestry showed that reduced 25(OH)D levels were associated with allcause mortality, cancer, and mortality from other causes [11].
Common genetic variants in ~ six loci have been identified to affect blood 25(OH)D concentrations in genomewide association studies (GWAS) performed in European whites [22], Hispanics [23], and Asian Indians [24]. Also, multiple GWAS on T2D have identified the association of common variants in over 100 genes in ethnically diverse cohorts [4]. However, whether genetically influenced reduced 25(OH)D concentrations can be causal or contributes to enhancing T2D risk, or conversely, whether gene variants that are now widely accepted to be robustly associated with increased risk for T2D could be involved in altering blood 25(OH)D levels is unclear and inconsistent.
Here, we have used genetic instrumental variable methods to obtain estimates of the causal association between circulating 25(OH)D concentrations and T2D, including the direction of causality, by performing two-directional reciprocal Mendelian randomization study [25]. The selection of index SNPs from six genetic loci as an instrumental variable was based on the bonafide gene-regions identified in GWAS and metanalysis studies for affecting T2D and blood 25(OH)D concentrations and their availability in all participant studies [22,24,[26][27][28]. These included two clusters of 3 index SNPs, each associated with 25(OH)D concentrations (GC rs2282679, DHCR7 rs12785878, CYP2R1 rs12794714), and T2D (IGF2BP2 rs1470579, TCF7L2 rs7903146, and KCNQ1 rs2237896), respectively. Even though multiple GWAS studies have pointed out the role of variants in genes involved in the synthesis, transport, and metabolism for influencing circulating 25(OH)D concentrations, the results of earlier published Mendelian randomization studies to determine the causal link between vitamin D insufficiency and T2D risk have been inconsistent and controversial. With few exceptions, these studies neither clearly support nor exclude the causal association of 25(OH)D levels with T2D [11,[29][30][31]. Moreover, these studies have examined mainly individuals from European populations, except one recent study included East Asians with Europeans [31]. Thus, in this study, we, for the first time,

Conclusion:
Our new data on this bidirectional Mendelianrandomization study suggests that genetically instrumented T2D risk does notcause changes in 25(OH)D levels. However, genetically regulated 25(OH)Ddeficiency due to vitamin D synthesis gene (DHCR7) may influence the risk ofT2D.
have included data from South Asian Indians with Europeans. People from the Indian sub-continent have 3 to sixfold higher T2D prevalence than Europeans and are a highly understudied population for vitamin D genetic studies. There is a paucity of data showing the relationship between vitamin D status and T2D in South Asian Indians [32].

Ethical statement
All participating studies were reviewed and approved by the respective Universities or Institutional Review Boards, including the primary institute of the University of Oklahoma Health Sciences Center's Institutional Review Board, as well as the Human Subject Protection Committees at the participating hospitals and institutes in India. All participants provided written informed consent for investigations.

Study design and phenotypic measurements
Our study with data from up to 59,890 participants from eight cohorts included: the Asian Indian Diabetic Heart Study/Sikh Diabetes Study (AIDHS/SDS); the Indian Migration Study (IMS); Twins UK; The British 1958 birth cohort (1958 BC); the Copenhagen City Heart Study (CCHS), the Copenhagen City General Population Study (CGPS), and the Copenhagen Ischemic Heart Disease Study (CIHDS); The UK Household Longitudinal Study (UKHLS); and The Prevention of Renal and Vascular End-Stage Disease (PREVEND) study. Details of study design and experimental work flow is presented in the Flow Chart (online Supplementary Section). Phenotypic measurements including T2D, 25(OH)D and other recruitment details are presented separately for each cohort (online Supplementary Section). Measures of blood 25(OH)D levels were not available on UKHLS (European) and MS India (Asian Indian) cohorts. UV index values were calculated for the city nearest the blood draw location and the month of recruitment using data from the National Weather Service Climate Prediction Center (http:// www. cpc. ncep. noaa. gov) and by averaging of the previous three months UV index as described previously [33]. The clinical characteristics of each study cohort stratified by T2D are summarized in Table 1.

SNP Genotyping and Quality Control (QC) and statistical analysis
The selection of SNPs as genetic instruments was based on their strong global association with respective exposure phenotypes of T2D and 25(OH)D reported in multiethnic GWAS and meta-analysis studies [22,[26][27][28]. Additional criteria for selecting instrumental variables were based on their availability in most of the participating studies. We have provided genotyping, sample, and SNP QC details separately by each cohort in the online Supplementary Section. Genome-wide genotype data for each cohort were checked for population structure using principal components, and outliers were removed before the analysis. To ensure that the Mendelian randomization assumptions were not violated, all SNPs selected as genetic instruments were strongly associated with the outcome (i.e. gene variants associated with T2D for increasing T2D outcome; and gene variants associated with vitamin D for affecting plasma 25(OH)D outcome). None of the SNPs used as Mendelian instruments for T2D and 25(OH)D were in linkage disequilibrium (r 2 = 0.001). None of the instrumental variables were associated with the outcome via exposure to other factors (BMI, age, and gender).
For 25(OH)D analysis, linear regression and additive genetic model were used with the natural log-transformed 25(OH)D level adjusted with age, gender, BMI, and study-specific covariates (e.g., T2D, UV index). Similarly, the three GWAS SNPs known to affect T2D risk were analyzed for their association with T2D using logistic regression and an additive genetic model. The T2D binary trait was the dependent variable with covariates age, gender, BMI. Each cohort analysis was limited to a single ancestral group (e.g. European ancestry only, South Asian only). If selected SNPs were unavailable, the best proxy SNP (SNP on the local haplotype having strong linkage disequilibrium (LD) with the reference SNP) was used using SNAP [34].
Upon completing the association analysis, the final list of SNPs was determined for gene score construction for each trait. The SNP, risk allele, number of risk alleles, and effect size (of each allele regarding the effect allele for log-transformed 25(OH)D levels) were required to construct a weighted gene score using the PLINK program [35]. Samples with missing 25(OH) D levels were excluded from the analysis (see Supplementary Section). A weighted gene score construction for T2D required the SNP, risk allele, number of risk alleles, and odds ratio (of each allele concerning the effect allele for T2D) for gene construction through PLINK as described previously [22]. Samples missing T2D status data were excluded, and SNPs that were directly genotyped are included to construct both gene scores. The additional parameter "-score-no-meanimputation" was used in PLINK to construct both gene scores. The "-score-no-mean-imputation" was used so that the missing allele was not accounted as the most common allele based on the sample allele frequencies while constructing gene scores. The constructed gene scores were each normalized to a quantitative variable at an appropriate scale.
An association analysis was performed on the constructed gene scores for each cohort. A linear regression and an additive genetic model were used with blood 25(OH)D levels as the dependent variable adjusted by gender, age, BMI, and T2D. Similarly, a logistic regression and an additive genetic model were used with T2D as the dependent variable adjusted by the gender, age, BMI. To combine the association analysis results of each trait by SNP for each cohort, a fixed effect inverse variance meta-analysis implemented in METAL [36] was performed. This required the minor allele frequency, effect size, standard error, and p-value for each SNP across the different ethnic cohorts. The combined results were stratified by their respective trait association; the three SNPs associated with 25(OH)D levels were assessed for their association with T2D and 25(OH)D. Similarly, the three significant SNPs for T2D were evaluated for their association with T2D and blood 25(OH)D. Each cohort provided their statistical data summaries for the SNPs, and meta-analysis was performed on these SNPs by trait using METAL and Forest Plot Viewer (http:// ntp. niehs. nih. gov/ ntp/ ohat/ fores tplot/), and PRISM (https:// www. graph pad. com/ scien tific-softw are/ prism/) was used to generate the forest plots. We did not observe heterogeneity of association of the genetic instruments of T2D and 25(OH) D with their respective phenotypes between studies; thus, fixed-effect models were used in metanalysis to drive odds ratios and confidence intervals. Finally, to assess the robustness of our conclusions, we applied conservative Bonferroni's correction and used a corrected significant threshold of 0.0063 (0.05/ 8 test variables including six instrument SNPs, T2D, and 25(OH) D concentration)

Results
Demographics and clinical characteristics of all study cohorts are summarized in Table 1. As expected, there were significant differences in the distribution of BMI and 25(OH)D among cases and controls in most participating cohorts. The overall distribution of blood 25(OH) D was consistent with previously published reports on these studies, and mean 25(OH)D levels were much lower across South Asian, compared to European cohorts. Data on 25(OH)D concentration was available for only a few individuals in the South Asian cohorts, except for the Punjabi Sikh studies (AIDHS/SDS).

Association of T2D susceptibility genes with T2D and blood 25(OH)D concentrations
The association of genetic variants with an increased risk for T2D was analyzed using three known and well-established T2D associated variants from previous studies. Meta-analysis of association of IGF2BP2 (rs1470579) with T2D (adjusted for age, gender, and BMI) showed an OR (95% CI) of 1  Table 2).

Genes involved in 25(OH)D synthesis, metabolism, and transport and their effects on 25(OH)D concentrations and T2D
Interestingly, not only was there a significant association, the direction of the minor (effect) allele of GC (rs2282679) and CYP2R1 (rs12794714) and DHCR7 (rs12785878) with blood 25(OH)D were largely consistent across chorts (Supplementary Table 3 Table 3). Next, we tested the association of the same variants with the risk of T2D. As summarized in Supplementary Table 4, neither variants for GC and CYP2R1 were associated with the risk of increasing T2D in any of the individual cohorts in the combined metaanalysis. The per-allele effect for increasing T2D risk yielded an OR (95% CI's) of 1.03 (0.97,1.09), P=0.207, and Z score = 0.16 for the GC variant (rs2282679), and 1.01 (0.97,1.05), P=0.453, Z score = 1.09 for the CYP2R1 (rs12794714). However, the per-allele effect of the DHCR7 (rs12785878) variant showed a significantly increased risk for T2D OR 1.05 95% CI (1.0, 1.11), P=0.004, and Z score =2.84 in the meta-analysis (Supplementary Table 4). We also adjusted the effects of UV index on modulating association of 25(OH)D SNPs on 25(OH)D concentrations (Supplementary Tables 5), and also for their association with T2D (Supplementary Tables 6). The overall outcome of association remained unchanged as the data on the UV index was only available in AIDHS/SDS.
Next, we combined the allelic effects of genes affecting blood 25(OH)D concentrations and gene variants associated with T2D risk by constructing a gene score. Table 4 shows a significant association between the gene score of 25(OH)D lowering alleles in GC (rs2282679), CYP2R1 (rs12794714), and DHCR7 (rs12785878) in a total of 41,136 individuals which yielded Beta ± SE (-0.021 ± 0.001),    Figure 2A). However, the same gene score of the 25(OH)D lowering allele was not associated with the increased risk for T2D OR, 1.002 (1.001-1.005), p= 0.12, and Z score = 1.54 (Table 5 and Supplementary Figure 2B).

Discussion
Here we performed a Mendelian randomization study to investigate whether genetically reduced vitamin D levels could be causally related to increasing T2D risk using data from up to 59,890 participants (5,862 cases and 54,028 controls). Essentially, our findings suggest that even though the candidate genetic variants (individually) or via composite gene score are associated with increased susceptibility to T2D, these are not involved in affecting blood 25(OH)D concentration. This means that the genetically instrumented T2D risk could be independent of the pathway (s) linked with vitamin D synthesis, transport, or early stages of metabolism. Additionally, our data showed that the genetically instrumented reduction   [11]. During exposure to sunlight, a precursor molecule, 7-dehydrocholesterol (DHC), is converted to cholesterol by the action of ultraviolet light (UV) on the skin through a thermal isomerization. Hydroxylation of pre-vitamin D3 occurs in the liver from the actions of CYPs (CYP2R1). The CYP27A1 drives the conversion of 25(OH)D to 1,25(OH)D in the kidney ( Supplementary  Fig. 3). While our present study was in progress, a study using the Mendelian randomization approach showed that a combined gene score of 25(OH)D genes involved in synthesis (DHCR7 rs12785878) and metabolism (CYP2R1 rs10741657) was associated with an increased risk for T2D in European and Chinese adults [31]. However, our results could not confirm the role of CYP2R1 variants in the risk for T2D. This discrepancy could be due to the use of a different CYP2R1 SNP (rs12794714) in the present study or due to the overlapping role of multiple cytochrome P450 (CYP) enzymes on the hydrolysis of 25(OH)D in hepatocytes [37].
As multiple CYPs in addition to CYP2R1 and CYP27A1 could evidently be involved in vitamin D hydroxylation in the liver [38]. Therefore, the role of other CYPs resulting in compensation for 25 hydroxylase activity cannot be ruled out [37]. Thus, in light of these findings, our negative findings on genetically reduced vitamin D by CYP2R1 are not surprising.
Presumably, the differences in linkage disequilibrium patterns among different populations, population-specific risk factors, and pleiotropy may have masked the potential cumulative effects of gene scores. It appears that perhaps due to similar limitations, the study of Lu et al. [31] also could not capture the impact of genetically regulated 25(OH)D concentration on T2D using a genetically instrumented gene score of four SNPs GC (rs2282679), CYP2R1 (rs10741657), Cyp27A1 (rs6013897) and DHCR7 (rs12785878) with T2D (p = 0.07) in meta-analysis despite using a sample size of up to 476,099 of Chinese and European adults. We undertook an additional meta-analysis combining the results of our analyses with those of Lu et al. [31]. These findings for the DHCR7 SNP showed that the genetically instrumented reduction in 25(OH)D of Beta ± SE (-4.3 nmol/L ± 0.3 nmol/L), p = 1.26 × 10 -98 , Z score -21.08 (Supplementary Table 7B) would significantly increase T2D risk from 5% (current study) to 7%; P = 5 × 10 -4 Z score = 3.51 in a total of 493,057 individuals (Supplementary Table 8B). However, we did not observe this effect in the meta-analysis of the GC variant with the results of Lu et al. [31] (Supplementary  Tables 7A and 8A). These data further confirm the possible causal effects associated with synthesis (DHCR) SNPs for increasing the risk of T2D.
The strengths of our study include well-characterized multiethnic cohorts to detect allelic association of vitamin D genes with T2D risk. We reduced the confounding effects of age, gender, and obesity across all studies. Principal components were used in each individual cohort to control for the potential of confounding by population stratification. Sample size-weighted Z score method was used to reduce inter-study variation in the diverse ancestries and differences in 25(OH)D assays between cohorts as described earlier [22]. Limitations of this study include insufficient data on 25(OH)D measures and T2D phenotypes in some cohorts; over-representation of European cohorts; lack of data on 25(OH)D measures in most South Asian studies that were available with GWAS; and data access was limited to summary statistics for each study cohorts. It is possible that the differences in serum 25(OH)D levels and T2D prevalence between South Asians and Europeans may influence the outcome of the association of biomarker with the genetic instrument. However, this would less likely be the case because; 1) the genetic instruments used for T2D and 25(OH)D were from bonafide candidate genes; 2) the application of sample-size weighted Z-score method would reduce inter-study variation, and 3) our metanalysis results do not show much heterogeneity in the SNP-phenotype association across all cohorts.
Additionally, seasonal differences in sun exposure (UV index) for effecting 25(OH)D concentration in most studies were not available except for our Punjabi Sikh cohort (AIDHS/SDS). Nonetheless, the geographic location of most participating European cohorts was nearby (Northern Europe), which might not have confounded the results due to latitude (to derive distance from the equator) variation (which all ranged from 54 0 to 55 0 N) (See Supplementary Table 10). Nevertheless, adjusting for the effects of the UV index in our Punjabi sample (AIDHS/ SDS) did not change the overall impact of genetically instrumented 25(OH)D on T2D susceptibility in AIDHS/ SDS and overall across studies.
It is possible to speculate that the lack of association with 25(OH)D metabolism and transport genes, but the association with synthesis gene (DHCR7) could suggest that any association with T2D may be through a UVdependent, vitamin D-independent effect as described [39]. The power of capturing the causal association would be improved by the expansion and the refinement of the appropriateness of the genetic instrument to determine causality. Moreover, because of some inherent limitations of Mendelian randomization design, including pleiotropy, buffering effects of environmental and genetic forces (redundancy due to more than one gene), and the effects of developmental compensation, large sample size may be required to logically confirm the causal association of 25(OH)D status with T2D.

Conclusions
For the first time, we have performed a bidirectional multiethnic Mendelian randomization study to determine the causal relationship between T2D and vitamin D concentrations, including data from the population of South Asians with Europeans. Even though our research has not found strong evidence of causation, it also does not rule out the possible contribution of genetic influence of vitamin D synthesis SNPs in increasing T2D risk. As the long-standing question of the role of vitamin D in influencing the risk of T2D remains unaddressed, our data stresses the need for the population-specific design for future observational and randomized clinical trials on this important and controversial topic of immense clinical importance.
In summary, our bidirectional Mendelian randomization study having data from South Asian Indians with Europeans suggests that genetically instrumented T2D risk may not be causing changes in 25(OH)D levels. However, we cannot entirely exclude the likelihood of DHCR7 genetic variants to influence the risk of T2D.

Supplementary Information
The online version contains supplementary material available at https:// doi. org/ 10. 1186/ s12937-021-00725-1.  (Lu et al, 2018) showing effects of genetically instrumented vitamin D candidate genes (GC and DHCR7) SNPs for their effects on T2D risk. Supplementary  Table 9. Information on samples with missing 25(OH)D levels from 8 different cohorts. Supplementary Table 10. Differences in latitude, 25(OH)D concentration and distribution of allele frequencies of vitamin D and T2D SNPs among South Asian and European cohorts.

Additional file 1: Supplementary
Additional file 2: Supplementary (Flow Chart). Experimental plan including the details of participating cohorts and candidate gene SNPs used as genetic instrument for performing a bidirectional Mendelian randomization study. Supplementary Figure 1 (A-B). Association analysis of composite gene score of T2D SNPs (IGF2BP2, TCF7L2 and KCNQ1) used as a genetic instrument for their joint effect on T2D and 25(OH)D concentrations (https:// www. graph pad. com/ scien tific-softw are/ prism/). Supplementary Figure 2 (A-B). Association analysis of composite gene score of vitamin D SNPs (GC, CYP2R1 and DHCR7) used as a genetic instrument for their effect on 25(OH)D concentrations and T2D (https:// www. graph pad. com/ scien tific-softw are/ prism/). Supplementary Figure 3. This figure illustrates pathways of synthesis, absorption, metabolism, and transportation of vitamin D. Genetic factors known to influence circulating levels of vitamin D are shown in the figure. DHCR7 expresses a reductase which uses nicotinamide adenine dinucleotide phosphate-oxidase to catalyzes the production of cholesterol to 7-dehydrocholesterol (7-DHC). The GC gene catalyzes the vitamin D binding protein formation. The enzyme (25-hydroxylase) responsible for the first hydroxylation step is encoded by CYP2R1. The enzyme in the kidneys responsible for the second hydroxylation is catalyzed by the CYP27B1 gene product (not used in this study). 1,25(OH)2D3 is the most effective and commonly measured vitamin D deficiency marker. Genetic instruments used in this study include gene pathways involved in vitamin D synthesis (DHCR7), metabolism (CYP2R1) and transport (GC) and genetic loci linked with increased T2D risk (IGF2BP2, TCF7L2, and KCNQ1). Note that in our results causality appears to be inferred by direct association of 25(OH)D variants (i.e. DHCR7) with 25(OH)D levels and the T2D phenotype. In contrast the T2D genes are associated with phenotype but not with the intermediate phenotype measured as 25(OH)D. UKHLS. The UK Household Longitudinal Study is led by the Institute for Social and Economic Research at the University of Essex. The UK Household Longitudinal Study is funded by the Economic and Social Research Council (ES/H029745/1). The survey was conducted by Nat Cen and the genome-wide scan data were analyzed and deposited by the Wellcome Trust Sanger Institute (WT098051). Information on how to access the data can be found on the Understanding Society website https:// www. under stand ingso ciety. ac. uk/.
IMS. This work was supported by the Wellcome Trust, UK (Grant no. 070797 and Grant no. 083541) which provided funding support for phenotyping and genotyping.