Serum bilirubin and kidney function: a Mendelian randomization study.

Background: Further investigation is needed to determine the causal effects of serum bilirubin on the risk of chronic kidney disease (CKD).
Methods: This study is a Mendelian randomization (MR) analysis. Among the well-known single-nucleotide polymorphisms (SNPs) related to serum bilirubin levels, rs4149056 in the SLCO1B1 gene was selected as the genetic instrument for single-variant MR analysis, as it was found to be less related to possible confounders than other SNPs. The association between genetic predisposition for bilirubin levels and estimated glomerular filtration rate (eGFR) or CKD was assessed in 337 129 individuals of white British ancestry from the UK Biobank cohort. Two-sample MR based on summary-level data was also performed. SNPs related to total or direct bilirubin levels were collected from a previous genome-wide association study and confounder-associated SNPs were discarded. The independent CKDGen meta-analysis data for CKD were employed as the outcome summary statistics.
Results: The alleles of rs4149056 associated with higher bilirubin levels were associated with better kidney function in the UK Biobank data. In the summary-level MR, both of the genetically predicted total bilirubin {per 5 µmol/L increase; odds ratio [OR] 0.931 [95% confidence interval (CI) 0.871-0.995]} and direct bilirubin [per 1 µmol/L increase; OR 0.910 (95% CI 0.834-0.993)] levels were significantly associated with a lower risk of CKD, supported by the causal estimates from various MR sensitivity analyses.
Conclusion: Genetic predisposition for higher serum bilirubin levels is associated with better kidney function. This result suggests that higher serum bilirubin levels may have causal protective effects against kidney function impairment.

INTRODUCTION

Chronic kidney disease (CKD) is a major comorbidity in modern medicine and is increasing in prevalence [1]. As kidney function is a pivotal factor that relates to various health outcomes, and because the socioeconomic burden of CKD itself is considerable, identifying factors that causally affect kidney function is an important health issue.

Bilirubin has been reported to be protective against kidney function impairment. Epidemiologic studies suggest that higher serum bilirubin levels, within physiologic ranges, are associated with better kidney prognosis [2-6] and similar findings have been reported for cardiovascular outcomes [7-12]. However, although previous studies have suggested these possibilities, it remains unclear whether bilirubin can serve as a therapeutic target for modulating kidney function [13-15]. This is mainly because no clinical trial modifying bilirubin is currently available and observational findings are inevitably affected by confounders and reverse causality.

Mendelian randomization (MR) is an approach that has been recently introduced in the medical field and has been used to identify important causal effects of various environmental and medical factors in complex diseases [16]. MR utilizes a genetic instrument that is minimally affected by confounders or reverse causation because the genotype is determined at birth. MR tests the association between genetic predisposition for exposure and health outcomes and since the randomization for the genotype is performed before birth, MR can provide causal estimates between complex exposures and outcomes [17]. In recent studies, MR has been introduced to reveal the causal factors related to kidney function traits [18-20].

In this study we aimed to reveal the causal effects of bilirubin on kidney function by performing MR in two population-scale databases. We hypothesized that higher genetically predicted bilirubin levels would be associated with better kidney function, suggesting the causal protective effect of serum bilirubin on the kidney.

MATERIALS AND METHODS

Ethical considerations

The study was performed in accordance with the Declaration of Helsinki. The study was approved by the institutional review boards of Seoul National University Hospital (no. E-1910-044-1067) and the UK Biobank consortium (application no. 53799). Since the study investigated materials from anonymous databases or summary-level data, the requirement for informed consent was waived by the institutional review boards.

Study setting

This study was an MR analysis in two population-scale databases (Figure 1). The study first used UK Biobank data to test the association between kidney function and genetic predisposition for elevated bilirubin via a single-nucleotide polymorphism (SNP) that is strongly associated with serum bilirubin levels. Replication was performed by a summary-level two-sample MR that implemented a set of SNPs previously reported from the UK Biobank data as the genetic instrument for MR of bilirubin levels by introducing independent outcome summary statistics.

FIGURE 1:

Study flow diagram. The study consisted of two parts: a single-variant MR analysis based on the individual-level data of the UK Biobank and a two-sample summary-level MR with the summary statistics from the CKDGen GWAS meta-analysis. As the UK Biobank and the CKDGen GWAS are independent, the two analyses were performed for the purpose of replication.

UK Biobank database

The UK Biobank is a prospective population-based cohort of >500 000 individuals 40–69 years of age gathered from 2006 to 2010 in the UK. We used the information to test the association between genetic variance in serum bilirubin level and kidney function, identify genetic variants associated with possible confounders and construct a genetic instrument consisting of a set of SNPs for summary-level MR. The details of the database have been published previously [21-23].

In the analysis, we included unrelated individuals of white British ancestry from the UK Biobank data with data passing the basic quality control filter. Those who were outliers in terms of heterozygosity or missing rate and those with sex chromosome aneuploidy were excluded. A total of 337 129 individuals were ultimately included in the genetic analysis. The details of the clinical information collection in the population are described in the Supplementary data, Methods.

Genetic instrument for the single-variant MR

We first tested the association between kidney function and a genetic variant that plays a decisive role in serum bilirubin level. Previous genome-wide association studies (GWASs) have consistently reported that two loci, uridine diphosphate-glucuronosyltransferase 1–1 (UGT1A1) [24-26] and solute carrier organic anion transporter family member 1B1 (SLCO1B1) [25, 27], are strongly associated with serum bilirubin levels. Among the SNPs in these genes, we selected rs887829 and rs4149056 as potential genetic instruments for serum bilirubin level, as in another recent MR study [28]. The variant rs887829 in the UGT1A1 gene has been reported to explain ˃30% of the variance in serum bilirubin level and is in nearly complete linkage disequilibrium with the genetic polymorphisms underlying Gilbert's syndrome. The SNP rs4149056 is in the SLCO1B1 gene that encodes the protein that transports bilirubin from the blood into the liver; this allele also explains a proportion of the variance of serum bilirubin level.

As MR requires meeting the assumption of ‘independence’ [16], namely, that the genetic instrument is not associated with any other confounder, we carefully investigated whether the variants could be associated with major confounders for kidney function in the UK Biobank data. The potential confounders included hypertension; diabetes mellitus; obesity; alanine aminotransferase level, which is known to be relatively specific for liver injury compared with other liver enzymes [29]; and serum albumin level. We performed logistic or linear regression with the allele status of the genetic variant for each possible confounder adjusted for age, sex, the first 10 principal components and the five possible confounders. If the significance of a regression reached a two-sided P-value <.05, the exposure variant was eliminated from consideration as a genetic instrument to be tested for its association with kidney function in the UK Biobank data.

In addition, as MR requires meeting the assumption of ‘relevance’ [16], we tested the association between the genetic instrument and serum bilirubin level in the UK Biobank data, estimating the explained variance by the multiple regression method and calculating the F statistic, which should be >10 to avoid weak instrument bias [30]. The genetic data analysis was performed using PLINK 2.0 (version alpha 2.3) and R (version 4.0.2; R Foundation for Statistical Computing, Vienna, Austria) [31].

Kidney function outcome for single-variant MR

The estimated glomerular filtration rate (eGFR) calculated by the Chronic Kidney Disease Epidemiology Collaboration (CKD-EPI) method in the UK Biobank data was used as the kidney function outcome [32, 33]. Since eGFR based on the measured creatinine values may be biased by diet or body composition, we examined the eGFR based on the serum cystatin C level. We also considered stage 3–5 CKD as an outcome, which included those with an eGFR <60 mL/min/1.73 m2 or a prevalent history of kidney replacement therapy, identified by self-reports or electronic admission records. The association between the genetic variant and eGFR or CKD was tested by linear or logistic regression analysis adjusted for age, sex and the first 10 principal components. In a sensitivity analysis, since eGFR values may have been altered from dialysis or transplantation, we reperformed the analysis in a population after additionally excluding 646 cases with confirmed kidney replacement therapy history.

Genetic instrument for the two-sample MR with summary-level data

To replicate the findings with an MR analysis based on summary-level data, we performed a two-sample MR analysis including a set of SNPs strongly associated with serum bilirubin levels. The approach of combining a number of SNPs as the genetic instrument has advantages over single-variant MR in terms of increased statistical power and ability to test for the presence of pleiotropy, thus the method was utilized for an independent outcome dataset for replication [34].

The genetic instrument used for the analysis was introduced in a previous large-scale GWAS for serum biomarkers in individuals of white British ancestry in the UK Biobank (Global Biobank Engine, https://biobankengine.stanford.edu/) [35]. The instrument has been utilized several times to genetically predict serum biomarkers [36-39]. The study identified genome-wide significant (P < 5 × 10–8) SNPs without linkage disequilibrium (R2 < 0.1) associated with total bilirubin or direct bilirubin levels.

To meet the independence assumption, the associations between the selected SNPs and five possible confounders, as above, were investigated by a GWAS with either linear or logistic regression, adjusted for age, sex and the first 10 principal components. As we aimed to robustly exclude SNPs that might be associated with the confounders, we applied a more stringent cutoff (P < .01) value to filter out confounder-associated SNPs than was used in previous studies [40, 41]. SNPs that showed a potential association with any of the confounders were excluded from the genetic instrument.

To meet the relevance assumption, we calculated the allele scores for serum bilirubin levels with the genetic variants after excluding the confounder-associated SNPs and then tested the significance of the association between the allele scores and the serum bilirubin level by linear regression analysis, adjusted for age, sex and the first 10 principal components. We also analyzed the explained variance by the multiple regression method with the calculation of F statistics.

In addition, we scaled the betas of the genetic instruments for total and direct bilirubin so that a unit of allele score reflected a 5 µmol/L (0.292 mg/dL) increase in total bilirubin and a 1 µmol/L (0.059 mg/dL) increase in direct bilirubin levels, respectively.

Summary statistics for CKD in the two-sample MR

The CKDGen consortium provides the largest database of GWAS meta-analysis results for stage 3–5 CKD (https://ckdgen.imbi.uni-freiburg.de) [42]. In the meta-analyzed data including 480 698 individuals of European ancestry, the prevalence of stage 3–5 CKD was ∼9%. Since the genetic instrument was developed from the individuals of white British ancestry in the UK Biobank data, we downloaded the summary statistics for CKD of the European ancestry individuals and utilized the data as the outcome statistics in our two-sample MR analysis.

Statistical method for the two-sample MR based on summary-level data

In the summary-level MR, any SNPs that did not overlap between the summary statistics or that were palindromic with intermediate allele frequencies were discarded [43]. The main method for the two-sample MR was the fixed-effects inverse variance weighted method. Sensitivity analyses were performed to calculate robust causal estimates independent of possible heterogeneity or pleiotropy. First, MR-Egger regression, which yields pleiotropy-robust causal estimates, was performed, with bootstrapped standard errors [44]. Second, the penalized weighted median mode method, which derives valid causal estimates even in conditions when invalid instruments are present, was implemented [45]. Finally, MR-pleiotropy residual sum and outlier (PRESSO), which detects and corrects the effects from outliers, yielding causal estimates that are robust to heterogeneity, was performed [46]. The two-sample MR analysis was performed by the TwoSampleMR package in R [47].

RESULTS

Baseline characteristics of the UK Biobank participants

The baseline characteristics of the 337 129 individuals of white British ancestry in the UK Biobank data included for the genetic analysis are described in Table 1. The median age was 58 years, with a 20.9% prevalence of hypertension and a 4.8% prevalence of diabetes mellitus. The interquartile ranges (IQRs) for laboratory values, including aspartate aminotransferase, alanine aminotransferase, albumin and direct or total bilirubin values, were all within the reference range. The prevalence of stage 3–5 CKD was 4.7%.

Table 1.

Baseline characteristics of the UK Biobank study population used for the genetic analysis

Characteristics	Total (n = 337 129)	Male (n = 156 106)	Female (n = 181 023)
Age (years), median (IQR)	58 (51–63)	59 (51–64)	58 (51–63)
Hypertension, n (%)	70 018 (20.9)	38 538 (24.9)	31 480 (17.5)
Systolic BP (mmHg), median (IQR)	136.5 (125.0–149.5)	139.5 (129.0–152.0)	133.5 (121.5–147.5)
Diastolic BP (mmHg), median (IQR)	82.0 (75.5–89.0)	84.0 (77.5–90.5)	80.0 (73.5–87.0)
Diabetes mellitus, n (%)	16 178 (4.8)	10 012 (6.4)	6166 (3.4)
Hemoglobin A1c (mmol/L), median (IQR)	35.1 (32.7–37.7)	35.2 (32.7–37.9)	35.1 (32.7–37.6)
Obesity (BMI ≥30 kg/m2), n (%)	81 022 (24.1)	39 328 (25.3)	41 694 (23.1)
BMI (kg/m2), median (IQR)	26.7 (24.1–29.8)	27.3 (25.0–30.0)	26.1 (23.4–29.6)
Laboratory values, median (IQR)
Aspartate aminotransferase (U/L)	20.2 (15.4–27.4)	23.8 (18.4–31.8)	17.5 (13.9–23.0)
Alanine aminotransferase (U/L)	24.4 (21.0–28.8)	26.1 (22.6–30.9)	23.0 (20.0–26.8)
Albumin (g/L)	45.2 (43.5–46.9)	45.5 (43.8–47.2)	45.0 (43.3–46.7)
Direct bilirubin (µmol/L)	1.6 (1.3–2.1)	1.8 (1.4–2.3)	1.5 (1.2–1.8)
Total bilirubin (µmol/L)	8.1 (6.4–10.4)	9.2 (7.4–11.7)	7.3 (5.9–9.2)
eGFR, cystatin C (mL/min/1.73 m2), median (IQR)	88.9 (77.1–101.3)	87.9 (76.7–99.7)	89.9 (77.5–102.6)
CKD stage 3–5, n (%)	15 004 (4.7)	7132 (4.8)	7872 (4.6)

BP, blood pressure; BMI, body mass index.

The genetic instrument for single-variant MR

In the analysis to test the independence assumption within the UK Biobank data, the allele status of rs887829 on UGT1A1 showed a certain association with the presence of diabetes mellitus and the presence of the T allele was associated with higher odds of diabetes (Supplementary data, Table S1). The rs4149056 SNP in SLCO1B1 did not show any significant association with any possible confounders. Thus rs4149056 was used as the genetic instrument to investigate the association between genetically predicted serum bilirubin level and eGFR.

When we tested the relevance assumption, rs4149056 was strongly (P < 2 × 10–16) associated with both total and direct bilirubin levels (for direct bilirubin: 0.4% of variance explained and F statistic 489.4; for total bilirubin: 0.5% of variance explained and F statistic 823.5). We confirmed that the C allele of the rs4149056 SNP was related to higher direct and total serum bilirubin levels (Supplementary data, Table S2).

Results of the single-variant MR

The presence of the C allele at the rs4149056 SNP was associated with higher eGFR values (Table 2). When we assessed stage 3–5 CKD as an outcome, the rs4149056 SNP was significantly associated with the risk of CKD, as the presence of the C allele, which was associated with higher bilirubin levels, was significantly associated with a lower risk of stage 3–5 CKD. The findings were similarly identified even if we excluded individuals with a history of kidney replacement therapy (Supplementary data, Table S3).

Table 2.

Genotypes of rs4149056 SNP and direct bilirubin, total bilirubin and eGFR values

		eGFR		CKD stage 3–5
Genotype (rs4149056)	eGFR (mL/min/1.73 m2), median (IQR)	Adjusted β (SE)	P-value	Adjusted OR (95% CI)	P-value
T/T	88.8 (70.0–101.2)	Reference		Reference
T/C	89.2 (77.3–101.5)	0.323 (0.063)	2.45 × 10–7	0.954 (0.918–0.992)	.018
C/C	90.0 (78.1–102.1)	1.092 (0.183)	2.39 × 10–9	0.882 (0.785–0.992)	.037

SE, standard error; OR, odds ratio; CI, confidence interval.

Adjusted for age, sex and the first 10 principal components.

The genetic instrument for summary-level MR

In the summary-level MR, among the 307 and 233 SNPs with dbSNP reference cluster (rs) identifiers that were related to total bilirubin and direct bilirubin levels, respectively, 212 and 126 SNPs were excluded from the genetic instrument in the summary-level MR as being potentially associated with at least one of the possible confounders (Supplementary data, Table S4). After additionally excluding SNPs for lack of overlap with the summary statistics of the CKDGen consortium and for being palindromic with intermediate allele frequencies, 79 and 66 SNPs remained eligible for inclusion in the genetic instrument for total and direct bilirubin, respectively (Supplementary data, Tables S5 and S6).

When we tested the relevance assumption by the allele scores calculated from the genetic instrument, allele scores for both total and direct bilirubin levels were strongly associated with phenotypic total and bilirubin levels, respectively (P < 2 × 10–16). The explained variance (total bilirubin 22.3%, direct bilirubin 14.3%) and F statistics (total bilirubin 9169, direct bilirubin 4564) also indicated that the genetic instruments were valid in regards to their strength of association with the phenotypes of interest.

Results of the summary-level MR

Genetic predispositions for both total and direct bilirubin were significantly associated with a lower risk of CKD according to the inverse variance weighted method (Figure 2 with OR and Table 3 with regressed betas and number of statistics). The MR-Egger pleiotropy test P-values (with a genetic instrument for total bilirubin: .175; for direct bilirubin: .358) and Cochran's Q statistics (with a genetic instrument for total bilirubin: .512; for direct bilirubin: .963) indicated that no significant pleiotropy or heterogeneity was likely to have been present in the causal estimates. Further, the causal estimates were supported by the sensitivity analysis results, as MR-Egger and penalized weighted median methods also indicated a significant association between genetically predicted bilirubin levels and the risk of CKD. We performed an MR-PRESSO analysis, but the global test for heterogeneity indicated no need to correct for heterogeneity (with a genetic instrument for total bilirubin: P = .571; for direct bilirubin: P = .974), so the outlier-corrected causal estimates remained the same as the raw results.

FIGURE 2:

Two-sample MR analysis for the causal estimates from the genetic predisposition for bilirubin on the risk in CKD. We scaled the betas of the genetic instruments for total and direct bilirubin so that a unit of allele score reflected a 5 µmol/L (0.292 mg/dL) increase in total bilirubin and a 1 µmol/L (0.059 mg/dL) increase in direct bilirubin levels, respectively. OR, odds ratio; CI, confidence interval.

Table 3.

The summary-level MR results for the causal estimates from serum bilirubin levels on risk of CKD

Exposure trait	MR method	Cochran's Q statistic for heterogeneity	MR-Egger pleiotropy test P-value	OR (95% CI)	P-value
Total bilirubin	Inverse variance weighted (fixed effects)	0.512	.175	0.931 (0.871–0.995)	.035
	MR Egger (bootstrap)			0.913 (0.849–0.982)	.004
	Penalized weighted median			0.919 (0.844–0.999)	.048
	MR-PRESSO[a]			0.931 (0.871–0.996)	.037
Direct bilirubin	Inverse variance weighted (fixed effects)	0.963	.358	0.910 (0.835–0.992)	.035
	MR Egger (bootstrap)			0.890 (0.809–0.980)	.008
	Penalized weighted median			0.894 (0.806–0.992)	.033
	MR-PRESSO[a]			0.910 (0.847–0.979)	.015

OR, odds ratio; CI, confidence interval.

The exposure betas were scaled such that a unit increase in allele score for total or direct bilirubin reflects a 5 µmol/L increase in total or a 1 µmol/L increase in direct bilirubin value.

MR-PRESSO global test to detect heterogeneity indicated the absence of significant heterogeneity (for total bilirubin: P = .571; for direct bilirubin: P = .974), therefore the causal estimates were the same as the raw analysis results.

DISCUSSION

In this MR study, we found that a genetic predisposition for higher serum bilirubin levels was significantly associated with higher eGFR and a lower risk of CKD. The results were repetitively identified in the individual-level data of the UK Biobank and also in the summary-level for the CKDGen GWAS. With our efforts to meet the necessary assumptions for an MR analysis, the study results support that serum bilirubin may have effects toward improved kidney function.

Beneficial effects of serum bilirubin on kidney and cardiovascular outcomes have been suggested by previous observational findings. The findings were first reported in individuals with Gilbert syndrome, who have a nonpathologic elevation of serum bilirubin and a reduced risk of adverse cardiovascular outcomes relative to the general population [7, 12]. Several cohort studies reported that higher serum bilirubin was associated with a lower risk of kidney function impairment in a wide range of ethnic populations [2–4, 9], which was confirmed by a systematized meta-analysis [6]. The possible association between higher serum bilirubin and lower risk of kidney function impairment was supported by experimental findings [2, 48]. Based on the above findings, bilirubin has been considered a potential therapeutic target for kidney and cardiovascular diseases.

However, due to the possibility of confounding or reverse causation effects, previous observational findings could not prove the benefits of serum bilirubin on human kidney function. This doubt was even enhanced by recent MR studies reporting the absence of a causal effect of serum bilirubin on cardiovascular diseases [13-15]. In addition, there have been no large-scale cohort studies assessing the causal effects of serum bilirubin levels on kidney function to date. MR analysis is a tool to estimate the causal effects of complex exposure on a health outcome. Through this study, we found that genetic predisposition for higher serum bilirubin is significantly associated with better kidney function parameters. Our study has strengths in that we performed a large-scale analysis and that the findings were consistent in both the single-variant MR and the summary-level MR, supporting that higher serum bilirubin may have a causal role in better kidney function.

MR requires that three assumptions are met in order to reveal the causal effects between complex exposures and diseases [16]. Through careful inspection of the genetic instrument to ensure the independence assumption was met, we excluded the SNPs that were possibly associated with potential confounders. In addition, in the summary-level MR, statistical tests indicated that no significant heterogeneity or pleiotropy biased the causal estimates. Furthermore, the relevance assumption was met and the genetic instruments utilized in this study were strongly associated with the phenotypic bilirubin levels. Although the remaining exclusion-restriction assumption cannot be formally tested, the weighted median method eases this assumption for up to half of the instrumented weight, again yielding significant causal estimates in the summary-level MR [45]. Thus our study made certain efforts to attain the key assumptions of MR and the findings suggest that serum bilirubin levels may be causally linked to kidney function.

There are several limitations that should be considered when interpreting our study results. First, as the study was based on a general population cohort, the study results do not indicate that a pathologic increase in serum bilirubin level would benefit kidney function. Second, currently there are few drugs that target human bilirubin levels, so there would be many obstacles to testing the actual benefits of interventions to increase serum bilirubin on kidney function. In addition, the MR result has limited usefulness to prove the benefits of relevant clinical intervention [49]. Third, the genetic analysis is not robust to detect nonlinear effects or to quantitatively estimate causal effects, thus, to what extent the bilirubin level is beneficial for kidney function cannot be answered by this study. Lastly, the study cohort comprises individuals of European ancestry, which limits the generalizability of these findings.

In conclusion, a genetic predisposition for higher serum bilirubin is significantly associated with better kidney function outcome. The MR findings support that higher serum bilirubin may be a protective factor for kidney function impairment. Further study is warranted to confirm the possible benefits of bilirubin modification on kidney function.

Click here for additional data file.