Circulating bilirubin, a natural antioxidant, is associated with decreased risk of type 2 diabetes (T2D), but the nature of the relationship remains unknown. We performed Mendelian randomization in a prospective cohort of 3,381 participants free of diabetes at baseline (age 28–75 years; women 52.6%). We used rs6742078 located in the uridine diphosphate–glucuronosyltransferase locus as an instrumental variable (IV) to study a potential causal effect of serum total bilirubin level on T2D risk. T2D developed in a total of 210 participants (6.2%) during a median follow-up period of 7.8 years. In adjusted analyses, rs6742078, which explained 19.5% of bilirubin variation, was strongly associated with total bilirubin (a 0.68-SD increase in bilirubin levels per T allele; P < 1 × 10−122) and was also associated with T2D risk (odds ratio [OR] 0.69 [95% CI 0.54–0.90]; P = 0.006). Per 1-SD increase in log-transformed bilirubin levels, we observed a 25% (OR 0.75 [95% CI 0.62–0.92]; P = 0.004) lower risk of T2D. In Mendelian randomization analysis, the causal risk reduction for T2D was estimated to be 42% (causal OR for IV estimation per 1-SD increase in log-transformed bilirubin 0.58 [95% CI 0.39–0.84]; P = 0.005), which was comparable to the observational estimate (Durbin-Wu-Hausman χ2 test, P for difference = 0.19). These novel results provide evidence that an elevated bilirubin level is causally associated with the risk of T2D and support its role as a protective determinant.
Introduction
Experimental and human studies show that elevated bilirubin levels are associated with decreased risk of type 2 diabetes (T2D) and diabetes-related outcomes (1–6). Bilirubin might antagonize oxidative stress by acting as an antioxidant and cytoprotectant, which may have beneficial effects in diseases related to oxidative stress (7–10). For example, we have previously shown that elevated bilirubin levels are associated with a reduction in mortality in the carriers of iron overload genotypes, which are phenotypically associated with increased oxidative stress (11). It is, however, unclear whether the association between bilirubin and diabetes is free of unobserved confounders. Apart from confounding, bilirubin concentrations could be changed as a consequence of T2D pathology (i.e., reverse causality). Therefore, causal inference of the effect of bilirubin on T2D (i.e., bilirubin being an intermediate phenotype in the etiology of T2D) is uncertain in conventional observational studies (12).
One approach to overcome these aforementioned limitations is to implement an instrumental variable (IV) that represents bilirubin levels and allows for assessment of the relation to T2D that is robust to confounding and reverse causality biases (13,14). Recent genome-wide association studies (GWAS) have identified that rs6742078, the single nucleotide polymorphism (SNP) that is mapped to hepatic uridine diphosphate-glucuronyltransferase (UGT1A1), explains 18% of the variation in serum total bilirubin levels (15–17). The UGT1A1 gene encodes hepatic uridine diphosphate-glucuronyltransferase, which is a major enzyme in bilirubin conjugation and in the regulation of bilirubin levels (5,15–17). When UGT1A1 is downregulated, bilirubin levels subsequently rise. Given the fact that the inheritance of either a T or a G allele of rs6742078 is randomly defined prior to disease onset (i.e., at the time of conception, based on the laws of Mendel), the variation of rs6742078 is robust to unobserved confounding and reverse causality. Therefore, rs6742078 is a suitable candidate to be used as an IV representing bilirubin variation between individuals. The approach of using genetic instruments (i.e., rs6742078) to test the causal association of a given intermediate exposure (i.e., bilirubin) with a disease outcome (i.e., T2D) is called Mendelian randomization (14). We hypothesized that implementing rs6742078 as an IV could provide evidence that the previously reported association between serum bilirubin level and T2D is causal or is not.
We aimed to investigate whether elevated total bilirubin levels are causally associated with a lower risk of new-onset T2D. By implementation of Mendelian randomization, we examined whether rs6742078(T) is associated with serum total bilirubin levels and T2D risk in a population-based cohort study.
Research Design and Methods
Study Population and Design
The Prevention of Renal and Vascular End-stage Disease (PREVEND) study is a Dutch population-based cohort for which 8,592 participants (age range 28–75 years) from the city of Groningen, the Netherlands, were recruited between 1997 and 1998. We examined a random sample of 4,016 individuals of the total cohort in which genotyping was performed. Details of the PREVEND study have been published previously (18). Among 4,016 individuals, 367 participants were excluded because of not passing quality control of DNA analyses, 124 participants were excluded because of missing follow-up data or information about medication, and 144 individuals were excluded because of prevalent diabetes at baseline, leaving 3,381 participants eligible for analyses (Fig. 1).
The PREVEND study has been approved by the medical ethics committee, University Medical Center Groningen, and conformed to the principles outlined in the Declaration of Helsinki. All participants gave written informed consent.
Baseline Study Measurements
At baseline and three rounds of screening from 1997 to 1998 (i.e., baseline examination) until 1 January 2007 (i.e., the third examination), study participants underwent two outpatient visits to obtain clinical characteristics and to collect two 24-h urine samples on 2 consecutive days (18). Furthermore, information on medication use was substantiated using pharmacy-based data from all community pharmacies in the city of Groningen (19). Blood samples were obtained after overnight fasting and were stored at −80°C until analyses without previous thawing and refreezing (Supplementary Table 1). Total bilirubin level was measured by a colorimetric assay (2,4-dichloroaniline reaction; MEGA; Merck, Darmstadt, Germany).
Definition of Outcome
During follow-up, participants were classified as new-onset T2D case patients if the fasting plasma glucose level was ≥7.0 mmol/L (126 mg/dL), the random sample plasma glucose concentration was ≥11.1 mmol/L (200 mg/dL), they reported a physician diagnosis of T2D, or they received insulin or oral hypoglycemic agents based on a central pharmacy registration. We included case patients in whom T2D was diagnosed from the first 3 months after the baseline screening visits (1997–1998) until January 2007 (20,21).
Genotyping and Quality Control
Genotyping was performed using Illumina HumanCytoSNP-12 arrays. Genotypes of 3,649 individuals were available preimputation (see Fig. 1 for quality control of samples). SNPs were called using Illumina Genome Studio software. SNPs were excluded with a minor allele frequency <0.01, a call rate <0.95, or a deviation from Hardy-Weinberg equilibrium (P < 1 × 10−5). Genome-wide genotype imputation was performed using Beagle version 3.1.0. We used the National Center for Biotechnology Information build 36 of Phase II HapMap CEU data (release 22) as the reference panel (22). Given the expected strong relationship of rs6742078 (imputation accuracy score >0.9) with bilirubin levels (16,17) and the fact that rs6742078 is in high linkage disequilibrium (LD) (r2= 0.88) with the functional UGT1A1*28 polymorphism (17), we decided that rs6742078 is a suitable IV for bilirubin in our study.
Statistical Analysis
We performed analyses in four steps. First, we examined the associations of bilirubin with participants’ characteristics, including age, sex, smoking status, family history of diabetes, and associated clinical variables. This was performed using linear regression. Next, we used logistic regression to examine the association of bilirubin with new-onset T2D (model A). This model estimates an observed effect of bilirubin on T2D (βbilirubin-T2D coefficient). The odds ratio (OR) and its corresponding 95% CIs for the risk of T2D was calculated per 1-SD increase in log-transformed bilirubin levels. Model A was adjusted for age, sex, family history of diabetes, BMI, hypertension, smoking status, and fasting glucose level (21). Second, we tested whether rs6742078 was associated with log-transformed bilirubin levels. We used linear regression under the assumption of an additive genetic model (model B), which is in agreement with previously published studies (17,23). Model B was adjusted for the same covariates as those included in model A (14). This analysis provided the observed effect estimate, which was expressed as the change in the number of SDs in log-transformed bilirubin per each copy increase in the number of T allele (range 0–2). To estimate the proportion of bilirubin variation explained by the IV of the rs6742078 genotypes, we calculated the difference between r2 values (i.e., explained variance) for model B without the genotype and for model B with the genotype. Third, logistic regression was used to evaluate the association between rs6742078 and T2D (model C). We calculated the OR for T2D per each copy increase in the rs6742078(T) allele. Model C was adjusted for the same covariates as those of model A. Fourth, we used the Wald-type method of IV estimation to test a potential causal relationship between elevated total bilirubin levels and decreased T2D risk. For the Wald-type method, the estimated risk is calculated using the following equation:
In this equation, ORgenotype-T2D is the OR resulting from logistic regression between rs6742078 and T2D, representing the risk of T2D per each copy of the rs6742078(T) allele (model C). The βgenotype-bilirubin coefficient is obtained from the regression of bilirubin on the SNP per each copy of the T allele of rs6742078 (model B). To validate the causal estimate derived from the Wald-type method, we used two alternative methods, called the two-stage least squares (2SLS) method and the multiplicative generalized method of moments (MGMM) (Table 3). The first stage of the 2SLS method is to examine the observational association between rs6742078 and bilirubin by means of linear regression (model 1), saving the predicted values and the residuals. In the second stage, the predicted values of bilirubin from model 1 are used as covariates, with T2D as the dependent variable in a logistic regression model, providing the βIV2SLS estimate (14). To obtain the MGMM βIV, the Stata command “ivpois” was used (24). The MGMM does not require distributional assumptions about bilirubin and thus can provide consistent estimation of causal OR (24).
To examine whether the association between bilirubin and T2D has “endogeneity,” we compared the IV estimate with the observed estimate from standard logistic regression analysis (i.e., βbilirubin-T2D coefficient) using the Durbin-Wu-Hausman χ2 test, which assesses whether there is a difference between these two estimates (14,24). We calculated the power of our study to detect a causal effect for bilirubin, given our sample size and bilirubin variation explained by rs6742078 (17,25).
To account for potential confounding bias, we tested whether rs6742078 is associated with common diabetes risk factors, other variables (Table 2), and 24-h urine albumin excretion (UAE). The association with UAE was tested because of the enrichment of the PREVEND cohort with microalbuminuric individuals at baseline (18).
To test for pleiotropy, we performed logistic regression with bilirubin as an independent variable and T2D as a dependent variable, and saved the residuals. Next, we tested the correlation between the genotype and the residuals that are independent of total bilirubin levels (26). To further dissect potential pleiotropy, we also calculated IV estimates from two additional bilirubin SNPs, namely, rs4149056 in SLCO1B1 and rs16928809 in SLC22A18 (17), in the PREVEND study. Additionally, we extracted effect estimates for rs6742078, rs4149056, and rs16928809 from publically available GWAS data on bilirubin, T2D (DIAbetes Genetics Replication And Meta-analysis [DIAGRAM]), or glycemic traits (Meta-Analyses of Glucose and Insulin-related traits Consortium [MAGIC]) (17,27–29) (Supplementary Table 2), and we calculated a summary statistics genetic risk score for bilirubin, as described previously (27). We then evaluated the presence of heterogeneity (derived from the likelihood-based method [27]) on the causal effects of these three SNPs on T2D, glucose, insulin, and HOMA of insulin resistance (HOMA-IR) (27–29) (Table 3). The presence of heterogeneity indicates potential pleiotropy (30).
We also performed a series of secondary analyses. First, we repeated regression models for bilirubin levels and rs6742078 by using a weighted method to compensate for the enrichment of the PREVEND cohort with microalbuminuric individuals (18). Second, we added the interaction terms of bilirubin × sex and bilirubin × age to model A or rs6742078 × sex and rs6742078 × age to model C, to evaluate potential sex or age differences in the association of bilirubin or rs6742078 with T2D, respectively. Third, we examined whether the effect estimates for bilirubin or rs6742078 were consistent across age strata of people <50 or ≥50 years old. Fourth, we accounted for the time between obtaining the sample for the measurement of bilirubin (and determining the genotype) and the development of T2D by using a Cox regression in which we applied the approximate time to the incidence of T2D as the time scale in our analysis. The time to the incidence of T2D cannot be ascertained exactly, because we screened for incident T2D periodically. We tested the proportional hazards assumption indicating that the effect estimates do not vary over time. Fifth, a potential correlation between the dropout rates and the genotype was tested. Sixth, we examined the association of rs6742078 with prevalent or a combination of both prevalent and incident T2D. A single imputation and predictive mean matching method was applied for missing data (21). For most baseline clinical variables, <1% of data was missing, whereas this was up to 8% for self-reported variables such as family history of diabetes. All the statistical analyses were carried out using IBM SPSS version 19.0, Stata/IC version 11.2 (StataCorp, College Station, TX), or R version 2.15.1 (Vienna, Austria) for Windows (http://cran.r-project.org/).
Results
Baseline characteristics of the population are shown in Table 1. During a median follow-up period of 7.7 years (Q1–Q3 7.4–8.0 years), T2D developed in a total of 210 individuals (6.2%) (from 1997 to 1998 until 2007).
Variables . | Total (n = 3,381) . | Nonconverters (n = 3,171) . | Incident T2D case patients (n = 210) . | P value . |
---|---|---|---|---|
Age (years) | 49.4 ± 12.4 | 49.0 ± 12.4 | 55.2 ± 9.9 | <0.001 |
Male sex | 1,710 (50.6) | 1,573 (49.6) | 137 (65.2) | <0.001 |
Current smoker | 1,198 (35.4) | 1,116 (35.2) | 82 (39.0) | 0.13 |
Ex-smoker | 1,240 (36.7) | 1,158 (36.5) | 82 (39.0) | |
Alcohol drinker | 2,573 (76.5) | 2,416 (76.6) | 157 (75.1) | 0.62 |
History of CVD | 134 (4.1) | 119 (3.9) | 15 (7.5) | 0.01 |
Family history of diabetes | 527 (15.6) | 468 (14.8) | 59 (28.1) | <0.001 |
BMI (kg/m2) | 26.0 ± 4.1 | 25.8 ± 4.0 | 29.3 ± 4.6 | <0.001 |
Waist circumference (cm) | 88.6 ± 13.1 | 87.9 ± 12.8 | 99.5 ± 11.7 | <0.001 |
Systolic blood pressure (mmHg) | 123.9 ± 19.0 | 123.3 ± 18.9 | 133.1 ± 18.6 | <0.001 |
Diastolic blood pressure (mmHg) | 71.8 ± 9.9 | 71.5 ± 9.9 | 76.1 ± 9.6 | <0.001 |
Hypertension | 936 (27.7) | 833 (26.3) | 103 (49.0) | <0.001 |
Fasting glucose (mmol/L) | 4.7 ± 0.6 | 4.7 ± 0.6 | 5.6 ± 0.8 | <0.001 |
HDL cholesterol (mmol/L) | 1.32 ± 0.40 | 1.34 ± 0.40 | 1.09 ± 0.29 | <0.001 |
Triglycerides (mmol/L) | 1.16 (0.84–1.68) | 1.13 (0.83–1.65) | 1.64 (1.17–2.45) | <0.001 |
Fasting insulin (mU/L) | 7.9 (5.5–11.9) | 7.7 (5.4–11.5) | 13.0 (8.5–21.4) | <0.001 |
Bilirubin (μmol/L) | 7.7 ± 3.6 | 7.7 ± 3.6 | 7.1 ± 3.1 | <0.001 |
AST (units/L) | 25.7 ± 9.5 | 25.6 ± 9.5 | 27.4 ± 8.2 | 0.009 |
ALT (units/L) | 23.8 ± 16.6 | 23.4 ± 16.6 | 29.6 ± 16.4 | <0.001 |
GGT (units/L) | 23.0 (16.0–38.0) | 23.0 (16.0–36.0) | 36.0 (24.7–51.7) | <0.001 |
Albumin (g/L) | 45.8 ± 2.8 | 45.8 ± 2.8 | 45.2 ± 3.2 | 0.003 |
ALP (units/L) | 65.3 ± 22.3 | 65.0 ± 22.5 | 70.0 ± 18.3 | 0.002 |
hs-CRP (mg/L) | 1.29 (0.58–2.95) | 1.24 (0.56–2.87) | 2.22 (1.19–4.65) | <0.001 |
UAE (mg/24 h) | 10.4 (6.6–20.3) | 10.2 (6.5–19.7) | 17.2 (8.8–36.5) | <0.001 |
Variables . | Total (n = 3,381) . | Nonconverters (n = 3,171) . | Incident T2D case patients (n = 210) . | P value . |
---|---|---|---|---|
Age (years) | 49.4 ± 12.4 | 49.0 ± 12.4 | 55.2 ± 9.9 | <0.001 |
Male sex | 1,710 (50.6) | 1,573 (49.6) | 137 (65.2) | <0.001 |
Current smoker | 1,198 (35.4) | 1,116 (35.2) | 82 (39.0) | 0.13 |
Ex-smoker | 1,240 (36.7) | 1,158 (36.5) | 82 (39.0) | |
Alcohol drinker | 2,573 (76.5) | 2,416 (76.6) | 157 (75.1) | 0.62 |
History of CVD | 134 (4.1) | 119 (3.9) | 15 (7.5) | 0.01 |
Family history of diabetes | 527 (15.6) | 468 (14.8) | 59 (28.1) | <0.001 |
BMI (kg/m2) | 26.0 ± 4.1 | 25.8 ± 4.0 | 29.3 ± 4.6 | <0.001 |
Waist circumference (cm) | 88.6 ± 13.1 | 87.9 ± 12.8 | 99.5 ± 11.7 | <0.001 |
Systolic blood pressure (mmHg) | 123.9 ± 19.0 | 123.3 ± 18.9 | 133.1 ± 18.6 | <0.001 |
Diastolic blood pressure (mmHg) | 71.8 ± 9.9 | 71.5 ± 9.9 | 76.1 ± 9.6 | <0.001 |
Hypertension | 936 (27.7) | 833 (26.3) | 103 (49.0) | <0.001 |
Fasting glucose (mmol/L) | 4.7 ± 0.6 | 4.7 ± 0.6 | 5.6 ± 0.8 | <0.001 |
HDL cholesterol (mmol/L) | 1.32 ± 0.40 | 1.34 ± 0.40 | 1.09 ± 0.29 | <0.001 |
Triglycerides (mmol/L) | 1.16 (0.84–1.68) | 1.13 (0.83–1.65) | 1.64 (1.17–2.45) | <0.001 |
Fasting insulin (mU/L) | 7.9 (5.5–11.9) | 7.7 (5.4–11.5) | 13.0 (8.5–21.4) | <0.001 |
Bilirubin (μmol/L) | 7.7 ± 3.6 | 7.7 ± 3.6 | 7.1 ± 3.1 | <0.001 |
AST (units/L) | 25.7 ± 9.5 | 25.6 ± 9.5 | 27.4 ± 8.2 | 0.009 |
ALT (units/L) | 23.8 ± 16.6 | 23.4 ± 16.6 | 29.6 ± 16.4 | <0.001 |
GGT (units/L) | 23.0 (16.0–38.0) | 23.0 (16.0–36.0) | 36.0 (24.7–51.7) | <0.001 |
Albumin (g/L) | 45.8 ± 2.8 | 45.8 ± 2.8 | 45.2 ± 3.2 | 0.003 |
ALP (units/L) | 65.3 ± 22.3 | 65.0 ± 22.5 | 70.0 ± 18.3 | 0.002 |
hs-CRP (mg/L) | 1.29 (0.58–2.95) | 1.24 (0.56–2.87) | 2.22 (1.19–4.65) | <0.001 |
UAE (mg/24 h) | 10.4 (6.6–20.3) | 10.2 (6.5–19.7) | 17.2 (8.8–36.5) | <0.001 |
Data are shown as the mean ± SD or median (Q1–Q3) for continuous variables, which were compared by independent t tests or Mann-Whitney U tests as appropriate. Categorical variables are shown as n (%) and were compared by χ2 test. ALP, alkaline phosphatase; ALT, alanine aminotransferase; AST, aspartate aminotransferase; GGT, γ-glutamyl transferase.
Association Between Bilirubin and T2D
We observed that total bilirubin levels were lower in individuals in whom T2D developed versus those in whom T2D did not develop (7.1 ± 3.1 vs. 7.7 ± 3.6 μmol/L, respectively; P < 0.001). Figure 2 depicts a graphical representation of the distribution of bilirubin with an incremental decrease in ORs for T2D over the range of bilirubin levels (2–35 μmol/L). The multivariable adjusted OR for T2D was 0.75 (95% CI 0.62–0.92), per 1-SD increase in log-transformed bilirubin levels (model A). The associations of bilirubin with baseline variables are shown in Table 2. When adding separately the variables, which were associated with bilirubin, to model A, corresponding ORs for T2D ranged from 0.78 (95% CI 0.67–0.92) to 0.83 (95% CI 0.69–0.94).
Variables . | Serum bilirubin . | rs6742078 . | ||
---|---|---|---|---|
Standardized β-coefficients (95% CI) . | P value . | Standardized β-coefficients (95% CI) . | P value . | |
Age | 0.012 (−0.021 to 0.046) | 0.46 | 0.013 (−0.021 to 0.046) | 0.46 |
Female sex | −0.215 (−0.245 to −0.181) | <0.001 | 0.001 (−0.033 to 0.034) | 0.97 |
Smoking | −0.015 (−0.050 to 0.017) | 0.35 | 0.003 (−0.030 to 0.037) | 0.87 |
Alcohol use | 0.049 (0.015–0.082) | 0.01 | 0.02 (−0.013 to 0.054) | 0.24 |
Family history of diabetes | −0.057 (−0.093 to −0.026) | 0.001 | 0.018 (−0.016 to 0.051) | 0.30 |
History of CVD | −0.026 (−0.060 to 0.007) | 0.13 | 0.026 (−0.003 to 0.064) | 0.14 |
BMI | −0.119 (−0.152 to −0.085) | <0.001 | 0.006 (−0.028 to 0.039) | 0.74 |
Waist circumference | −0.025 (−0.059 to 0.008) | 0.14 | −0.002 (−0.036 to 0.032) | 0.91 |
Systolic blood pressure | 0.039 (0.00–0.072) | 0.02 | −0.007 (−0.040 to 0.027) | 0.70 |
Diastolic blood pressure | 0.040 (0.003–0.070) | 0.03 | −0.003 (−0.036 to 0.031) | 0.88 |
Hypertension | 0.013 (−0.020 to 0.046) | 0.44 | −0.017 (−0.051 to 0.016) | 0.32 |
Fasting glucose | −0.021 (−0.055 to 0.013) | 0.22 | 0.025 (−0.009 to 0.058) | 0.15 |
HDL cholesterol | 0.025 (−0.009 to 0.058) | 0.15 | −0.006 (−0.039 to 0.027) | 0.73 |
Triglycerides | −0.116 (−0.148 to −0.082) | <0.001 | 0.001 (−0.033 to 0.034) | 0.97 |
Fasting insulin | −0.106 (−0.140 to −0.072) | <0.001 | −0.002 (−0.036 to 0.032) | 0.87 |
AST | 0.145 (0.112–0.178) | <0.001 | −0.003 (−0.036 to 0.030) | 0.85 |
ALT | 0.082 (0.049–0.115) | <0.001 | 0.021 (−0.013 to 0.054) | 0.22 |
GGT | 0.030 (−0.007 to 0.060) | 0.12 | 0.010 (−0.023 to 0.043) | 0.55 |
Albumin | 0.247 (0.212–0.282) | <0.001 | −0.014 (−0.052 to 0.023) | 0.44 |
ALP | −0.022 (−0.054 to 0.013) | 0.24 | 0.022 (−0.015 to 0.060) | 0.23 |
hs-CRP | −0.115 (−0.147 to −0.080) | <0.001 | −0.010 (−0.044 to 0.023) | 0.55 |
UAE | −0.031 (−0.066 to 0.002) | 0.07 | 0.006 (−0.027 to 0.040) | 0.72 |
Variables . | Serum bilirubin . | rs6742078 . | ||
---|---|---|---|---|
Standardized β-coefficients (95% CI) . | P value . | Standardized β-coefficients (95% CI) . | P value . | |
Age | 0.012 (−0.021 to 0.046) | 0.46 | 0.013 (−0.021 to 0.046) | 0.46 |
Female sex | −0.215 (−0.245 to −0.181) | <0.001 | 0.001 (−0.033 to 0.034) | 0.97 |
Smoking | −0.015 (−0.050 to 0.017) | 0.35 | 0.003 (−0.030 to 0.037) | 0.87 |
Alcohol use | 0.049 (0.015–0.082) | 0.01 | 0.02 (−0.013 to 0.054) | 0.24 |
Family history of diabetes | −0.057 (−0.093 to −0.026) | 0.001 | 0.018 (−0.016 to 0.051) | 0.30 |
History of CVD | −0.026 (−0.060 to 0.007) | 0.13 | 0.026 (−0.003 to 0.064) | 0.14 |
BMI | −0.119 (−0.152 to −0.085) | <0.001 | 0.006 (−0.028 to 0.039) | 0.74 |
Waist circumference | −0.025 (−0.059 to 0.008) | 0.14 | −0.002 (−0.036 to 0.032) | 0.91 |
Systolic blood pressure | 0.039 (0.00–0.072) | 0.02 | −0.007 (−0.040 to 0.027) | 0.70 |
Diastolic blood pressure | 0.040 (0.003–0.070) | 0.03 | −0.003 (−0.036 to 0.031) | 0.88 |
Hypertension | 0.013 (−0.020 to 0.046) | 0.44 | −0.017 (−0.051 to 0.016) | 0.32 |
Fasting glucose | −0.021 (−0.055 to 0.013) | 0.22 | 0.025 (−0.009 to 0.058) | 0.15 |
HDL cholesterol | 0.025 (−0.009 to 0.058) | 0.15 | −0.006 (−0.039 to 0.027) | 0.73 |
Triglycerides | −0.116 (−0.148 to −0.082) | <0.001 | 0.001 (−0.033 to 0.034) | 0.97 |
Fasting insulin | −0.106 (−0.140 to −0.072) | <0.001 | −0.002 (−0.036 to 0.032) | 0.87 |
AST | 0.145 (0.112–0.178) | <0.001 | −0.003 (−0.036 to 0.030) | 0.85 |
ALT | 0.082 (0.049–0.115) | <0.001 | 0.021 (−0.013 to 0.054) | 0.22 |
GGT | 0.030 (−0.007 to 0.060) | 0.12 | 0.010 (−0.023 to 0.043) | 0.55 |
Albumin | 0.247 (0.212–0.282) | <0.001 | −0.014 (−0.052 to 0.023) | 0.44 |
ALP | −0.022 (−0.054 to 0.013) | 0.24 | 0.022 (−0.015 to 0.060) | 0.23 |
hs-CRP | −0.115 (−0.147 to −0.080) | <0.001 | −0.010 (−0.044 to 0.023) | 0.55 |
UAE | −0.031 (−0.066 to 0.002) | 0.07 | 0.006 (−0.027 to 0.040) | 0.72 |
ALP, alkaline phosphatase; ALT, alanine aminotransferase; AST, aspartate aminotransferase; GGT, γ-glutamyl transferase; hs-CRP, high-sensitivity C-reactive protein.
Association Between rs6742078 and Bilirubin
We examined whether rs6742078, which has been previously discovered as the top-ranked SNP in a GWAS of total bilirubin (16,17), was associated with bilirubin levels in our cohort. In model B, the rs6742078*T allele was strongly associated with bilirubin levels (a 0.68-SD increase in log-transformed bilirubin levels per allele, P < 1 × 10−122). The proportion of bilirubin variation explained by rs6742078 was 19.5% (the difference between r2 = 0.070 for model B without the genotype and r2 = 0.265 for model B with the genotype). Using linear regression analysis, the rs6742078*T allele was associated with a 2.607 μmol/L higher total bilirubin level (Fig. 3, left panel).
Association Between rs6742078 and T2D
We observed that a one-copy increase in the number of the rs6742078*T allele was associated with T2D risk with an OR of 0.69 (95% CI 0.54–0.90; P = 0.006), which corresponds to a 31% decreased risk of T2D (Fig. 3, right panel). Adjustment for bilirubin level attenuated the association between rs6742078 and T2D risk (OR 0.80 [95% CI 0.60–1.06]; P = 0.13).
Causal Estimates for the Effect of Bilirubin on T2D risk
To estimate the causal effect of bilirubin on T2D, we used the Wald-type method, which assesses the ratio of the association of rs6742078 with T2D and that of rs6742078 with bilirubin levels (Fig. 4). Based on this, we calculated that the estimated causal effect of total bilirubin (per 1-SD increase in log-transformed bilirubin levels) corresponds to a 42% lower risk of T2D (causal ORIVestimation 0.58 [95% CI 0.39–0.84]; P = 0.005). Given a sample size of 3,381 individuals (including 210 incident cases of T2D) and 19.5% of bilirubin variation explained by rs6742078 in our study, we had 80% power at a significance level of 0.05 to estimate an OR of 0.64 per 1-SD increase in log-transformed bilirubin levels. We observed that the direction of the IV and observed estimates were similar with a slightly larger effect size for the IV estimate (Pendogeneity = 0.19). The use of two alternative methods for IV estimation provided evidence of causal effect, which was also comparable to the observed OR for T2D (Pendogeneity = 0.18 or 0.22; Table 3).
IV . | Alleles (effect/other) . | Mendelian randomization analysis . | Heterogeneity P value . | ||
---|---|---|---|---|---|
Method . | Causal estimate (95% CI) . | P value . | |||
Incidence of T2D, PREVEND study per 1-SD log | |||||
rs6742078 | T/G | Wald-type | 0.58 (0.39–0.84)* | 0.005 | |
2SLS | 0.56 (0.38–0.81)* | 0.002 | |||
MGMM | 0.50 (0.27–0.92)* | 0.027 | |||
rs4149056 | C/T | MGMM | 0.58 (0.05–6.08)* | 0.65 | |
rs16928809 | A/G | MGMM | 0.17 (0.02–2.16)* | 0.17 | |
Genetic risk score | MGMM | 0.34 (0.15–0.77)* | 0.009 | ||
T2D, DIAGRAM study (28)† per 1-unit log (μmol/L) | |||||
rs6742078 | T/G | Summary statistics | 1.14 (0.96–1.35)* | 0.13 | |
rs4149056 | C/T | Summary statistics | 0.67 (0.28–1.62)* | 0.38 | |
rs16928809 | A/G | Summary statistics | 0.44 (0.15–1.35)* | 0.15 | |
Genetic risk score | Summary statistics | 1.09 (0.93–1.29)* | 0.28 | 0.14 | |
Glucose, MAGIC study (29) per 1-unit log (μmol/L) | |||||
rs6742078 | T/G | Summary statistics | −0.04 (−0.07 to −0.01)‡ | 0.02 | |
rs4149056 | C/T | Summary statistics | 0.01 (−0.19 to 0.20)‡ | 0.95 | |
rs16928809 | A/G | Summary statistics | −0.05 (−0.28 to 0.18)‡ | 0.67 | |
Genetic risk score | Summary statistics | −0.04 (−0.07 to −0.01)‡ | 0.02 | 0.89 | |
Insulin, MAGIC study (29) per 1-unit log (μmol/L) | |||||
rs6742078 | T/G | Summary statistics | −0.04 (−0.07 to −0.00)‡ | 0.03 | |
rs4149056 | C/T | Summary statistics | 0.12 (−0.08 to 0.32)‡ | 0.25 | |
rs16928809 | A/G | Summary statistics | −0.04 (−0.27 to 0.20)‡ | 0.76 | |
Genetic risk score | Summary statistics | −0.03 (−0.07 to −0.00)‡ | 0.05 | 0.33 | |
HOMA-IR, MAGIC study (29) per 1-unit log (μmol/L) | |||||
rs6742078 | T/G | Summary statistics | −0.04 (−0.08 to −0.00)‡ | 0.03 | |
rs4149056 | C/T | Summary statistics | 0.08 (−0.13 to 0.28)‡ | 0.47 | |
rs16928809 | A/G | Summary statistics | −0.06 (−0.31 to 0.19)‡ | 0.63 | |
Genetic risk score | Summary statistics | −0.04 (−0.07 to −0.00)‡ | 0.04 | 0.55 |
IV . | Alleles (effect/other) . | Mendelian randomization analysis . | Heterogeneity P value . | ||
---|---|---|---|---|---|
Method . | Causal estimate (95% CI) . | P value . | |||
Incidence of T2D, PREVEND study per 1-SD log | |||||
rs6742078 | T/G | Wald-type | 0.58 (0.39–0.84)* | 0.005 | |
2SLS | 0.56 (0.38–0.81)* | 0.002 | |||
MGMM | 0.50 (0.27–0.92)* | 0.027 | |||
rs4149056 | C/T | MGMM | 0.58 (0.05–6.08)* | 0.65 | |
rs16928809 | A/G | MGMM | 0.17 (0.02–2.16)* | 0.17 | |
Genetic risk score | MGMM | 0.34 (0.15–0.77)* | 0.009 | ||
T2D, DIAGRAM study (28)† per 1-unit log (μmol/L) | |||||
rs6742078 | T/G | Summary statistics | 1.14 (0.96–1.35)* | 0.13 | |
rs4149056 | C/T | Summary statistics | 0.67 (0.28–1.62)* | 0.38 | |
rs16928809 | A/G | Summary statistics | 0.44 (0.15–1.35)* | 0.15 | |
Genetic risk score | Summary statistics | 1.09 (0.93–1.29)* | 0.28 | 0.14 | |
Glucose, MAGIC study (29) per 1-unit log (μmol/L) | |||||
rs6742078 | T/G | Summary statistics | −0.04 (−0.07 to −0.01)‡ | 0.02 | |
rs4149056 | C/T | Summary statistics | 0.01 (−0.19 to 0.20)‡ | 0.95 | |
rs16928809 | A/G | Summary statistics | −0.05 (−0.28 to 0.18)‡ | 0.67 | |
Genetic risk score | Summary statistics | −0.04 (−0.07 to −0.01)‡ | 0.02 | 0.89 | |
Insulin, MAGIC study (29) per 1-unit log (μmol/L) | |||||
rs6742078 | T/G | Summary statistics | −0.04 (−0.07 to −0.00)‡ | 0.03 | |
rs4149056 | C/T | Summary statistics | 0.12 (−0.08 to 0.32)‡ | 0.25 | |
rs16928809 | A/G | Summary statistics | −0.04 (−0.27 to 0.20)‡ | 0.76 | |
Genetic risk score | Summary statistics | −0.03 (−0.07 to −0.00)‡ | 0.05 | 0.33 | |
HOMA-IR, MAGIC study (29) per 1-unit log (μmol/L) | |||||
rs6742078 | T/G | Summary statistics | −0.04 (−0.08 to −0.00)‡ | 0.03 | |
rs4149056 | C/T | Summary statistics | 0.08 (−0.13 to 0.28)‡ | 0.47 | |
rs16928809 | A/G | Summary statistics | −0.06 (−0.31 to 0.19)‡ | 0.63 | |
Genetic risk score | Summary statistics | −0.04 (−0.07 to −0.00)‡ | 0.04 | 0.55 |
The bilirubin genetic risk score consisted of rs6742078 plus two other bilirubin SNPs in SLCO1B1, rs4149056 and in SLC22A18, rs16928809 (17,27–29).
*The causal ORs for T2D were estimated using the Wald-type method, the 2SLS method, the MGMM, and the summary statistics (27,30). The causal OR for T2D was calculated as OR = exp(βIV).
†The DIAGRAM GWAS includes both prevalent and incident cases of T2D.
‡The causal β-coefficient (95% CI) for glucose (mmol/L), insulin, and HOMA-IR was estimated for each SNP separately and the genetic risk score, which is based on the combination of these three SNPs.
Confounding and Pleiotropy
We found no evidence that rs6742078 is associated with potential confounders like various demographic, lifestyle, and clinical variables (Table 2). Thus, we assumed that the association between the IV and T2D was mainly free of confounding. To test for pleiotropy, we assessed the correlation between rs6742078 and the residuals from the logistic regression for T2D (as dependent variable) and bilirubin (as independent variable). There was no evidence for pleiotropy (P for correlation = 0.11). We consistently observed no heterogeneity for the causal effects of bilirubin on T2D risk or glycemic traits when we used a summary statistics genetic score including rs6742078, rs4149056, and rs16928809 (Table 3). Details of the association of individual SNPs with bilirubin, T2D, and glycemic traits are given in Supplementary Table 2. In the PREVEND study, while rs6742078 explained 19.5% of bilirubin variation, rs4149056 and rs16928809 explained only 0.18% and 0.15% of the variance, respectively. This limits the use of these SNPs as a single IV. For rs4149056 and rs16928809, we estimated the causal ORs of 0.58 (95% CI 0.05–6.08) and 0.17 (95% CI 0.02–2.16) in the PREVEND study. These estimates lack power for the detection of associations as these two SNPs explained only a very small amount of bilirubin variation and therefore are to be considered weak IVs.
Secondary Analyses
In subsequent analyses, we repeated the analyses in models A and C, when weighted for individuals with a mean urine albumin level of >10 mg/L. The results of these complex design analyses were similar to those of the main analyses, with an OR of 0.79 (95% CI 0.66–0.95; P = 0.01) per 1-SD increase in log-transformed bilirubin, and an OR of 0.59 (95% CI 0.47–0.97; P < 0.001) for the rs6742078*T allele. Second, we found no evidence for interactions of sex and age with bilirubin and rs6742078 (P for sex interaction = 0.49 and P for age interaction = 0.38 with bilirubin; P for sex interaction = 0.96 and P for age interaction = 0.25 with rs6742078) for the risk of T2D. Third, we observed that both observational and IV estimates were age consistent, with ORs of 0.69 (95% CI 0.51–1.04) and 0.45 (95% CI 0.21–0.94) in people <50 years of age, and 0.86 (95% CI 0.70–1.04) and 0.55 (95% CI 0.34–0.90) in people ≥50 years of age, respectively. Fourth, the assumption of proportional hazards was met for the Cox regression model, indicating there is no evidence of time-varying associations of bilirubin and rs6742078 with T2D (P = 0.38 and P = 0.54, respectively). Fifth, we observed no evidence of an association between the dropout rates (82.9% and 71.3% of participants, respectively, who underwent the second and the third round) and the genotype (P for correlation = 0.15). Sixth, we found an OR of 0.93 (95% CI 0.69–1.26) per each T allele of rs6742078 for prevalent T2D; and when we performed a combined analysis of both prevalent and incident T2D, our findings remained materially unchanged (OR 0.77 [95% CI 0.64–0.93]). In further analysis, we tested for another UGT1A1 SNP, rs887829 (in perfect LD with rs6742078; r2 = 1.0), which has been shown to be strongly associated with bilirubin levels in GWAS data (31). In our cohort, the rs887829*T allele was associated with a decreased risk of T2D (OR 0.70 [95% CI 0.54–0.89]), which was comparable to that of the rs6742078*T allele. Finally, we observed no evidence of systematic differences between the participants who were included in the main analysis compared with the rest of the cohort participants without diabetes at baseline (data not shown).
Discussion
In this prospective cohort study, we investigated whether we could provide evidence that bilirubin is causally associated with the risk of new-onset T2D. We demonstrated that the rs6742078*T allele in the UGT1A1 locus is strongly associated with elevated levels of endogenous bilirubin and is also associated with new-onset T2D. Several studies (1–6) have previously demonstrated associations between bilirubin and T2D in animal models and in humans. However, confounding and reverse causality might have contributed to the implication of causal relationships in observational studies. To our knowledge, this is the first study to report a potential causal association between elevated total bilirubin levels and decreased T2D risk that is likely free of potential confounders using Mendelian randomization.
Accumulating evidence shows that bilirubin has powerful antioxidant properties (7–10). Bilirubin may compensate for the oxidative stress that might be an important factor in the pathophysiology of diabetes (32–36). Oxidative stress is consistently related to increased plasma glucose levels and may contribute to the development of diabetes (37), and to microvascular and macrovascular complications of diabetes (38,39). Bilirubin has been shown to prevent oxidative stress in a number of diseases, including atherosclerosis, cancer, and diabetic nephropathy (40,41). In mice, biliverdin, a precursor of bilirubin, protects against the deterioration of glucose tolerance (3), while increased levels of bilirubin reduced streptozotocin-induced pancreatic β-cell damage through attenuating oxidative stress (42) or increasing insulin sensitivity in mice (43). In line with experimental studies, previous observational studies (2,44) have shown that elevated bilirubin level is associated with a 26–31% decreased risk of diabetes in humans.
In previous studies (5,15–17), variations in the UGT1A1 and SLCO1B1 loci seemed to explain part of interindividual differences in bilirubin levels. Johnson et al. (17) found that rs6742078, very close to the UGT1A1*28 TATA box polymorphism, accounted for ∼18% of the variation in total bilirubin levels (P < 5.0 × 10−324). Less powerful, but still present, was the association with rs4149056, a functional SNP of the SLCO1B1 locus, which accounted for 0.6% of bilirubin variation (17). In our cohort study, we observed that 19.5% of bilirubin variation was explained by the UGT1A1 rs6742078 SNP. It is likely that the association of rs6742078 with bilirubin can be explained to a large extent by the functional UGT1A1*28 genetic variant, which has been found to be in high LD (r2 = 0.88) with rs6742078 (17). Therefore, the contribution of genetic variation in UGT1A1 to the risk of T2D could be greater than estimated from the associations observed for rs6742078. The genetically determined component of protective effects may also vary globally as the TA repeat in the UGT1A1 promoter has other, less common functional variants (e.g., the UGT1A1*6 variant, which is more often observed in Asians with Japanese or Chinese ancestry) (45).
From phenotypic perspective, the Framingham Heart Study (5) has found that the UGT1A1*28 variant was associated with both bilirubin level and cardiovascular disease (CVD). However, the relationships of the UGT1A1 SNPs with CVD, glycemic traits, and diabetes have been inconsistent, possibly due to differences in design and populations across studies (23). In the current study, we found that the rs6742078 SNP was associated with T2D risk. We tested whether rs6742078 was associated with glycemic traits that were recently evaluated in meta-analysis of 21 GWA studies (29). rs6742078 was associated with fasting glucose level (β = −0.009, P = 0.02), insulin level (β = −0.009, P = 0.03), and HOMA-IR (β = −0.009, P = 0.03) based on a priori selected single-test association analysis in published GWAS data (29). In a recent GWAS meta-analysis (28) (12,171 case subjects and 56,862 control subjects), no evidence of an association between rs6742078 and T2D has been reported (P = 0.1). The lack of association between rs6742078 and T2D might be due to the heterogeneity in the T2D phenotype (28,46), as T2D cases have been ascertained using different sources (28). Of note, we identified a similar pattern of absence of an association in the GWAS meta-analysis, while another prospective Mendelian randomization study (47) reported the presence of such an association. In the latter study, the presence of evidence for the association between two SNPs in the gene encoding sex hormone-binding globulin and T2D could not be confirmed in the case-control GWAS meta-analysis. Case-control studies typically included both prevalent and incident cases, which differ with respect to the disease stage, diagnostic criteria, and duration of diabetes (46). Also, these discrepancies may be partly explained by differences in data quality control and preanalysis preparations (28,47). A prospective cohort design and its accompanying statistical methods will reduce potential biases such as incidence-prevalence bias and nonresponse bias (46). Although we estimated the risk of future T2D, there still might be a chance of a type I error in our analysis. The current results, together with the associations of rs6742078 with glycemic traits in the meta-analysis of GWAS, support the idea that elevated bilirubin levels may lead to decreased risk of insulin resistance and T2D. Our finding remains to be confirmed for the potential role of bilirubin in the trajectories of glycemic traits and the development of T2D by prospective cohort studies.
In our study, we used a strong IV, rs6742078, which satisfied the requirements for Mendelian randomization (48), to examine the potential causal relationship between bilirubin and T2D in a large population-based study. Using this approach, we directly investigated a relationship between genotype and phenotype independent of potential confounders. Other bilirubin SNPs, rs4149056 and rs16928809 (17), may not be selected as a single IV, given the fact that the strength of these instruments is limited (25) because of the small amount of bilirubin variation that is explained by these two SNPs (0.18% and 0.15% in the PREVEND study, respectively).
Our findings may have implications for the prevention or prediction of T2D and its complications. For example, Gilbert syndrome, which is characterized by moderate hyperbilirubinemia, exists in 5–10% of the general population (5,16,49). It has been shown that diabetic patients with Gilbert syndrome have reduced markers of oxidative stress and decreased risk of diabetic nephropathy and CVD (5,49). One may speculate that the modulation of bilirubin levels (e.g., using medical interventions that inhibit UGT1A1) might have protective effects against the risk of T2D (45). A search for the sources of increase in bilirubin levels or its antioxidative capacity (e.g., effective therapeutics with minimal adverse effects) merits further investigations.
One major concern remaining in the context of Mendelian randomization is that we cannot completely control for the possibility of reverse causality, unobserved confounding, and pleiotropy (14,50). The functional UGT1A1*28 variant (in high LD with rs6742078), which is associated with a 10–35% reduction in UGT1A1 activity, is mainly involved in glucuronidation of bilirubin, but also of other compounds (e.g., estrogen) and drugs (e.g., irinotecan) (45,51). Yet, the effects of the genotype on phenotypic traits other than bilirubin and their consequences on disease outcomes (e.g., cancer) remain inconclusive. An extensive knowledge of gene function and associated biological processes is needed to better understand the extent to which such genetic variation in the human genome contributes to lifelong protection against T2D. Our cohort is predominantly comprised of white adults of Dutch descent, and it is therefore unknown whether the findings of our study can be generalized to other populations (18). The PREVEND cohort was enriched for individuals with microalbuminuria at baseline. However, in a secondary analysis, we used a weighted method to compensate for this, and this did not affect our results.
In summary, endogenous bilirubin and the most highly associated SNP, rs6742078 in the UGT1A1 locus, are associated with the risk of T2D. We observed consistent effect estimates derived from alternative genetic instruments, rs4149056 and rs16928809. However, these two SNPs explain far less of bilirubin variation, and therefore the estimates derived from these individual instruments should be interpreted cautiously. Our findings provide evidence that lifelong genetically elevated bilirubin level are likely to be causally protective against the development of T2D. Further studies are warranted to validate this finding and to elucidate the exact underlying mechanisms of action.
Article Information
Acknowledgments. The authors thank Professor Dr. L.T.W. de Jong-van den Berg and Dr. S.T. Visser from the Department of Social Pharmacy, Pharmacoepidemiology, and Pharmacotherapy, Groningen University Institute for Drug Exploration, University of Groningen, University Medical Center Groningen, for providing the data on pharmacy-registered use of insulin or oral hypoglycemic agents.
Funding. This work was supported by the Netherlands Heart Foundation, Dutch Diabetes Research Foundation, Dutch Kidney Foundation, the Netherlands Organization for Scientific Research project (NWO), and the Medical Research Council UK (grant no. MC_U106179471). A.A. is supported by a Rubicon grant from the NWO (project no. 825.13.004).
Duality of Interest. No potential conflicts of interest relevant to this article were reported.
None of the study sponsors had a role in the study design; data collection, analysis, and interpretation; report writing; or the decision to submit the report for publication.
Author Contributions. A.A., P.E.D., and S.J.L.B. helped to conceive and design the study, analyze the data, write the first draft of the manuscript, and agreed with the manuscript results and conclusions. E.C., R.T.G., R.O.B.G., H.L.H., P.v.d.H., and G.N. contributed to the writing of the manuscript, and agreed with the manuscript results and conclusions. R.P.S. helped to conceive and design the study, contributed to the writing of the manuscript, and agreed with the manuscript results and conclusions. B.Z.A. helped to conceive and design the study, analyzed the data, contributed to the writing of the manuscript, and agreed with the manuscript results and conclusions. A.A., B.Z.A., and S.J.L.B. are the guarantors of this work and, as such, had full access to all the data in the study and take responsibility for the integrity of the data and the accuracy of the data analysis.