Introduction

Circulating total bilirubin has been consistently shown to be inversely and independently associated with adverse cardiometabolic outcomes such as cardiovascular disease (CVD), hypertension and type 2 diabetes [1,2,3]. Though a causal association has been demonstrated for total bilirubin and type 2 diabetes [4], there is no strong evidence for a causal association between total bilirubin levels and CVD [5, 6]. Nonalcoholic fatty liver disease (NAFLD), emerging as the most common cause of chronic liver disease in the developed world [7], is a cardiometabolic condition which is characterized by hepatic steatosis with varying degrees of necroinflammation and fibrosis [7]. In the absence of the reference standard—liver biopsy [8], the diagnosis of NAFLD is commonly based on (1) imaging techniques [i.e., ultrasonography, computed tomography (CT) scan, or magnetic resonance imaging (MRI)] confirming the presence of fat infiltration of the liver, (2) exclusion of other liver diseases of other aetiology and (3) in the absence of substantial alcohol intake [9]. Due to the high costs associated with these imaging techniques and their unsuitability for use in large-scale population-based studies, a number of biomarker-based algorithms have been developed to aid the diagnosis of NAFLD. The fatty liver index (FLI) [10] and the hepatic steatosis index (HSI) [11] are based on easily accessible variables and have been reported to have good diagnostic accuracies for NAFLD [10,11,12].

Given the existence of a strong link between total bilirubin and adverse cardiometabolic outcomes, which has been attributed to the antioxidant [13, 14], anti-inflammatory [15] and antiatherogenic properties of bilirubin; [16] there have been suggestions that circulating total bilirubin might also be associated with a reduced risk of NAFLD. Indeed, a number of observational studies have reported on the associations between bilirubin and NAFLD over recent years. However, these were based on cross-sectional or case–control study designs limited by lack of temporality, paediatric populations, selected patients with pre-existing disease or were unable to demonstrate associations specifically between total circulating bilirubin and NAFLD risk [17,18,19,20,21,22]. Therefore uncertainty remains regarding the nature and magnitude of the prospective association between total bilirubin and NAFLD. With the ongoing debate on the potential value of using circulating total bilirubin to prevent and treat cardiometabolic outcomes [23,24,25], it will be clinically useful if circulating total bilirubin is shown to contribute to the development of NAFLD. In this context, we aimed to quantify the nature and magnitude of the prospective association between total bilirubin and the risk of NAFLD (as estimated by these two indices—FLI and HSI) in a general population sample who were free from pre-existing disease (NAFLD, CVD, and malignancy) at baseline. Since observational epidemiological studies are beset by several biases such as residual confounding, reverse causation, and regression dilution [26,27,28,29], we utilized a Mendelian randomisation (MR) to assess if there is a causal relevance to the association using the well-known rs6742078 variant in the UGT1A1 gene [6, 30].

Materials and methods

Study population

We conducted this study according to STROBE (STrengthening the Reporting of OBservational studies in Epidemiology) guidelines for reporting observational studies in epidemiology (“Appendix 1”) [31]. Participants for the current analysis were part of the ongoing Prevention of Renal and Vascular End-stage Disease (PREVEND) study, a large-scale general population-based observational cohort study which began in 1997 in the Netherlands. The PREVEND study was designed to investigate the natural course of urinary albumin excretion and its relationship to renal and CVD progression. The study design and recruitment procedures have been described in detail in several previous reports [2, 32]. Briefly, 8592 inhabitants aged 28–75 years living in the city of Groningen, the Netherlands were recruited into the PREVEND study with baseline measurements performed between 1997 and 1998. For the present analysis, we excluded participants (1) with a prevalent history of CVD, renal disease, malignancy, or NAFLD (for the analysis of FLI outcome, participants with prevalent NAFLD as measured by FLI were excluded and vice versa for HSI) and (2) with excessive alcohol use (defined as four or more drinks per day), which left a cohort of 3824 participants (based on FLI) with non-missing information on total bilirubin levels, relevant covariates, and outcomes. Of these participants, 1610 individuals had complete phenotypic and genotypic data. The PREVEND study complies with the Declaration of Helsinki and was approved by the medical ethics committee of the University Medical Center Groningen. All participants provided written informed consent for voluntary participation.

Risk factor assessment

For the assessment of baseline data on sociographic characteristics, physical measurements, medical history, and use of medication, participants completed two outpatient visits. Further information on medication use was complemented using data from all community pharmacies in the city of Groningen and this covers complete information on drug use in 95% of PREVEND participants. Venous blood was obtained from participants after an overnight fast and 15 min of rest. Plasma samples were prepared by centrifugation at 4 °C. Blood lipids (total cholesterol, high-density lipoprotein cholesterol (HDL-C), and triglycerides (TGs)), high sensitivity C-reactive protein (hsCRP), serum creatinine, and serum cystatin C were measured using standard laboratory protocols previously described [33,34,35]. Total bilirubin was measured using a colorimetric assay (2,4-dicholoraniline reaction; Merck MEGA, Darmstadt, Germany), with the detection limit being 1.0 µmol/L. The inter-assay coefficients of variation were 3.8% and 2.9% in the lower normal and higher normal range respectively. The mean of two 24-h urine collections was used to estimate urinary albumin excretion (UAE) and its concentration determined by nephelometry (BNII; Dade Behring Diagnostic, Marburg, Germany). Serum liver aminotransferase (alanine aminotransferase, ALT and aspartate aminotransferase, AST) activities were measured using the standardized kinetic method with pyridoxal phosphate activation (Roche Modular P; Roche Diagnostics, Mannheim, Germany). Serum gamma-glutamyltransferase (GGT) activity was measured by an enzymatic colorimetric method (Roche Modular P; Roche Diagnostics, Mannheim, Germany). Estimated glomerular filtration rate (eGFR), was calculated using the Chronic Kidney Disease Epidemiology Collaboration (CKD-EPI) combined creatinine-cystatin C equation [36]. Body mass index (BMI) was estimated by dividing weight measured in kilograms by the square of height in meters.

NAFLD ascertainment

The FLI was calculated based on the report by Bedogni et al. [10] using the following formula:

$${\text{FLI}} = [{\text{e}}^{{(0.953 \times \ln \left( {\text{TG}} \right) + 0.139 \times {\text{BMI}} + 0.718 \times \ln \left( {\text{GGT}} \right) + 0.053 \times {\text{WC}} - 15.745)}} ] / [1 + {\text{e}}^{{(0.953 \times \ln \left( {\text{TG}} \right) + 0.139 \times {\text{BMI}} + 0.718 \times \ln \left( {\text{GGT}} \right) + 0.053 \times {\text{WC}} - 15.745)}} ]\times 100$$

The FLI ranges from 0 to 100, with FLI < 30 ruling out (sensitivity = 87%) and FLI ≥ 60 ruling in fatty liver disease (specificity = 86%) with a good diagnostic accuracy of 0.84 (95% CI 0.81–0.87).

The HSI was estimated using the following formula as reported by Lee et al.: [11]

$$\begin{aligned} {\text{HSI}} & = 8 \times {\text{ALT}}/{\text{AST ratio}} + {\text{BMI }}\left( { + 2,{\text{if diabetes}}; + 2,{\text{if female}}} \right); \\ & \quad {\text{with HSI values}} < 30{\text{ and}} > 36{\text{ ruling out and ruling in fatty liver respectively}}. \\ \end{aligned}$$

Genotyping

Genotyping was performed using Illumina HumanCytoSNP-12 arrays. SNPs were called using Illumina Genome Studio software. SNPs were excluded with minor allele frequency < 0.01, call rate < 0.95, or deviation from Hardy–Weinberg equilibrium (p < 1×10 − 5). The rs6742078 SNP was a suitable instrumental variable for the present analyses, given its robust specificity for circulating total bilirubin levels (explaining up to 45% of the variation in circulating bilirubin levels [37]) and its use in previous studies to assess the causal relevance of total bilirubin to several disease outcomes [5, 6, 30].

Statistical analyses

Observational analyses

Skewed variables (total bilirubin, ALT, AST, TGs, hsCRP, and UAE) were natural log-transformed to approximate normal distributions. Descriptive analyses were used to summarize baseline characteristics of participants overall and according to NAFLD outcomes. Partial correlation coefficients adjusted for age and sex were calculated to assess cross-sectional associations of total bilirubin levels with risk markers for NAFLD. Logistic regression was used to examine the association of total bilirubin with new-onset NAFLD as measured by FLI or HSI (observed effect of bilirubin on NAFLD). The odds ratio (OR) with 95% confidence intervals (CIs) for the risk of NAFLD was calculated per 1 standard deviation, SD higher loge total bilirubin levels. In subsidiary analyses, total bilirubin was modeled as tertiles defined according to its baseline distribution. The SD of baseline loge total bilirubin level was 0.42 (equivalent to 1.5-fold higher circulating total bilirubin level, as e0.42 = 1.52). Odds ratios were progressively adjusted for in four models: (Model 1) age and sex; (Model 2) plus smoking status, systolic blood pressure (SBP), total cholesterol, and HDL-C; (Model 3) plus alcohol consumption, glucose, eGFR, and UAE; and (Model 4) plus hsCRP. Confounders were selected based on their known associations with NAFLD and observed associations with total bilirubin using the available data [38] and evidence from previous research [1, 2]. We used interaction tests to assess statistical evidence of effect modification by sex on the association. To minimise bias due to potential reverse causation, we carried out sensitivity analyses that excluded participants with a history of diabetes at baseline, participants on regular statin medication or participants with potential Gilbert’s syndrome (defined as defined as total bilirubin > 34.2 µmol/L, AST < 80 IU/L, ALT < 80 IU/L, and GGT < 80 IU/L) [1].

Mendelian randomization

To assess for pleiotropy, we tested whether rs6742078 was associated with relevant risk markers which might confound the relationship between total bilirubin and NAFLD. We assessed the association between rs6742078 and loge total bilirubin using linear regression under the assumption of an additive genetic model, with adjustment for covariates used in the observational analysis between total bilirubin and NAFLD. This analysis provided the number of SDs above loge total bilirubin per each copy increase in number of the T allele. In order to estimate the causal effect of bilirubin on NAFLD, we used the two-stage least squares (2SLS) method. The first stage of the 2SLS method is to examine the association between rs6742078 and total bilirubin by means of linear regression (model 1) and then saving the predicted values and the residuals. In the second stage, the predicted values of total bilirubin from model 1 are used as covariates, with NAFLD as a dependent variable in a logistic regression, which provides the MR estimate. Estimates of power for the MR analysis on NAFLD employed an online power calculator ([http:/cnsgenomics.com/shiny/mRnd/]). We used the genetic sample size and case/control ratios together with the proportion of variance of total bilirubin explained by rs6742078. Using PREVEND data, we had 20% power to detect a causal association of bilirubin on NAFLD, similar in magnitude to the fully adjusted observational OR of 0.82. All statistical analyses were conducted using Stata version 15 (Stata Corp, College Station, Texas, USA).

Results

Baseline characteristics and correlates of total bilirubin

Baseline descriptive characteristics of study participants as well as cross-sectional correlates of total bilirubin are shown in Table 1. The overall mean age of participants at study entry was 47 (SD 12) years and 40.6% were women. The mean (SD) of loge total bilirubin level was 1.96 (0.42) µmol/l. There were weak and inverse correlations of loge total bilirubin levels with physical measures (BMI and waist circumference), lipids (total cholesterol and TGs), and glucose. There were weak positive correlations with HDL-C and the liver aminotransferases. There was an inverse correlation with loge hsCRP (r = − 0.24). Baseline total bilirubin levels were higher by 28% in men compared with women. The levels were lower by 12% in the combined group of current and former smokers compared with non-current smokers (Table 2). “Appendices 2–3” show baseline characteristics according to the development of NAFLD. Except for history of diabetes and smoking status, there were significant differences in baseline clinically relevant subgroups and levels of risk markers between participants who did and did not develop NAFLD (for both indices) during follow-up.

Table 1 Baseline characteristics and cross-sectional correlates of total bilirubin
Full size table
Table 2 Association of baseline total bilirubin with incident NAFLD as measured by FLI
Full size table

Total bilirubin levels and risk of incident NAFLD

During a mean (SD) follow-up of 4.2 (0.4) years, 434 and 452 cases of NAFLD as estimated by FLI and HSI respectively, were recorded. The associations between total bilirubin and NAFLD as estimated by the FLI are reported in Table 2. In age and sex adjusted analysis, the OR for NAFLD (as estimated by FLI) per 1 SD change in loge total bilirubin was 0.72 (95% CI 0.64–0.80; p < 0.001), which was minimally attenuated to 0.77 (95% CI 0.69–0.87; p < 0.001) after further adjustment for several established and emerging risk factors. Following additional adjustment for hsCRP, the OR was 0.82 (95% CI 0.73–0.92; p < 0.001). In analyses that compared the top versus bottom tertiles of total bilirubin, the inverse associations between total bilirubin and NAFLD were evident (Table 2).

Odds ratios for the associations between total bilirubin and NAFLD as estimated by the HSI are reported in Table 3. The age- and sex-adjusted OR for NAFLD per 1 SD change in loge total bilirubin was 0.79 (95% CI 0.72–0.88; p < 0.001), which remained consistent 0.84 (95% CI 0.75–0.93; p = 0.001) following further adjustment for several established and emerging risk factors. Following additional adjustment for hsCRP, the OR was 0.87 (95% CI 0.78–0.97; p = 0.012). In analyses that compared the top versus bottom tertiles of total bilirubin, the inverse associations remained except for evidence of attenuation to the null on further adjustment for hsCRP (Table 2).

Table 3 Association of baseline total bilirubin with incident NAFLD as measured by HSI
Full size table

The association between total bilirubin and NAFLD was not statistically significantly modified by sex for both measures of NAFLD (p > 0.05); whereas the association between total bilirubin and NAFLD (as estimated by FLI) was comparable in males and females, the age-adjusted inverse association between total bilirubin and NAFLD (as estimated by HSI) in males was attenuated to null on further adjustment for established risk factors (Table 4). The ORs for all associations remained similar in sensitivity analyses that involved exclusion of participants with prevalent diabetes, participants on cholesterol lowering medication, or participants with potential Gilbert’s syndrome (“Appendices 4–5”).

Table 4 Association of baseline total bilirubin with incident NAFLD in males and females
Full size table

Mendelian randomization findings

There was strong evidence for an association between rs6742078 SNP and total bilirubin (0.70 SD increase in total bilirubin levels per T allele; SE = 0.03, p = 1.35 × 10−80) and the SNP explained 20% of total variance in total bilirubin levels. There was no evidence for associations between rs6742078 and confounders included in the observational analyses, except for a weak association with sex (“Appendix 6”). The causal OR for FLI was 0.98 (95% CI 0.69–1.38; p = 0.900) per 1 SD genetically elevated total bilirubin level and that for HSI was 1.14 (95% CI 0.81–1.59; p = 0.451) (Table 5).

Table 5 Causal estimates for NAFLD (as measured by FLI and HSI) using Mendelian randomisation analysis
Full size table

Comments

Key findings

Using a large-scale population-based study of predominantly Caucasian men and women without pre-existing disease including NAFLD at baseline, we have demonstrated that total bilirubin is inversely associated with future risk of NAFLD as estimated by two well-known biomarker indices, the FLI and HSI respectively. The associations were independent of several established risk factors and other potential confounders and remained robust in several sensitivity analyses. The inverse association between total bilirubin and NAFLD was not significantly modified by sex. However, given the rather relatively small numbers (low event rates) in males and females, larger-scale studies in both genders are needed to further evaluate a potential role of sex on the impact of bilirubin on new onset NAFLD. Furthermore, using a MR analysis, we investigated whether the observed association between total bilirubin and NAFLD was devoid of confounding and/or reverse causation. Despite strong observational evidence for associations between bilirubin and NAFLD, these results were not supported by MR analyses—there was no evidence for a causal relationship between genetically elevated bilirubin and NAFLD. While the rs6742078 SNP was strongly associated with levels of circulating bilirubin and the F statistic (> 400) indicated good instrument strength, the observed wide 95% CIs for the MR results are indicative of low power. For MR to be valid, several assumptions need to be met: [39] the genetic instrument (1) must be associated with the exposure of interest, (2) must not be associated with any confounders and (3) must not directly influence the outcome, except through the exposure of interest. When exploring the association between rs6742078 and confounding variables included in the observational analyses, no evidence for associations was found, except for weak evidence of an association between rs6742078 and sex. While this could be a chance finding, given the multiple tests that were performed, it is also possible that this might represent a true causal relationship. This might reflect the sex differences in hepatic uridine diphosphate glucuronyltransferase (UGT1A1) activity [40], the enzyme that contributes substantially to bilirubin glucuronidation and enhances bilirubin elimination. Consistent with our findings, previous studies have demonstrated higher levels of circulating bilirubin in men compared with women [19]. On the contrary, a number of studies have reported findings which suggest that the discordance in circulating bilirubin levels between males and females may not be attributable to differences in the UGT1A1 polymorphism [41,42,43]. Therefore, these findings deserve follow-up in larger cohort studies.

Comparison with previous studies

A limited number of epidemiological observational studies have reported on the associations of bilirubin with the risk of NAFLD. In an analysis of over 17,000 participants, Kwak et al. demonstrated an inverse association between total bilirubin and NAFLD diagnosed on the basis of ultrasonographic findings and an alcohol consumption of less than 20 g/day [19]. However, the main limitation of this study was its cross-sectional design, which precluded the ability to assess the temporal relationship between bilirubin and the risk of NAFLD. Findings from two prospective cohort studies have demonstrated serum direct bilirubin to be associated with reduced risk of NAFLD, but found no evidence of associations for total or indirect bilirubin [17, 20]. To our knowledge, this is the first study to demonstrate a long-term prospective association between total bilirubin and the risk of NAFLD in a general predominantly Caucasian population. We are also the first to investigate whether a potential causal association exists between elevated total bilirubin levels and decreased NAFLD risk in a general adult population. Furthermore, our study employed NAFLD outcomes diagnosed on the basis of biomarker-based indices, whereas previous relevant studies have utilized ultrasonographic findings [17,18,19,20].

Potential explanations for findings

Bilirubin, iron and carbon monoxide are degradation products of heme catabolism via heme oxygenase-1 (HO-1) and they have been reported to regulate important functions in cells [44]. Increasing evidence suggests that bilirubin has cytoprotective properties [45], antioxidant actions [13, 14] and anti-inflammatory effects via its anti-complement actions [15, 46]. Since oxidative stress mechanisms and inflammation play a major role in the pathophysiology of NAFLD [47, 48], if there is a protective effect of bilirubin on NAFLD risk, this may be via the potent antioxidant and anti-inflammatory effects of bilirubin. Whether total bilirubin is a direct cause of NAFLD or just a risk marker of underlying NAFLD, is not certain at the moment. Our MR findings did not provide evidence for a causal association between circulating total bilirubin and NAFLD. It may be possible that the observational association may be driven by biases such as unmeasured confounding and/or reverse causation. The current evidence suggests that total bilirubin may just be a risk marker of NAFLD. However, given that the MR estimates were under-powered and a causal effect of bilirubin on NAFLD cannot be completely ruled out, further investigation in a larger sample is necessary.

Implications of findings

The current findings further contribute to the accumulating body of evidence on the putative protective role of bilirubin on adverse cardiometabolic outcomes such as hypertension, diabetes, metabolic syndrome and CVD [1,2,3,4]. Indeed, the utility of using circulating bilirubin to prevent and treat cardiometabolic disease has been extensively debated [23,24,25]. It is reported that patients with Gilbert syndrome (characterised by moderate hyperbilirubinaemia) have reduced levels of markers of oxidative stress and inflammation and have decreased risk of vascular complications; [49, 50] hence, it is possible that the sustained hyperbilirubinaemia prevents vascular complications via inhibition of oxidative stress and inflammation. Circulating bilirubin is a simple, standardised, cheap, and easily measured biomarker, hence its potential value in reducing the risk of adverse cardiometabolic outcomes warrants further evaluation.

Strengths and limitations

In addition to the novelty, several strengths of this work deserve mention. These include employing a sample that was representative of the general population, exclusion of participants with pre-existing disease which minimised biases due to reverse causation, accounting for several key clinical characteristics, and conducting several sensitivity analyses to confirm the robustness of the findings. In our MR analysis, we used the rs6742078 variant to examine the potential causal relationship, which has been shown to be strongly associated with substantial increases in levels of bilirubin and explained substantial proportion of variation in bilirubin levels. The limitations of the current work include the inability to generalise the findings to other populations as the sample predominantly comprised of white adults of Dutch descent, absence of data on repeat measurements of bilirubin to correct for regression dilution bias, and inability to replicate or verify the findings in an independent cohort due to lack of data. One might say that the association between rs6742078 and sex violated one of the conditions for the MR analysis, hence impacting on the causal estimates. However, this is unlikely for the following reasons: (1) the association was weak and (2) it is very likely this was a chance finding (due to the multiple statistical tests) given that the rs6742078 variant has consistently been demonstrated to be very specific for bilirubin and is not associated with various demographics (including sex), lifestyle and clinical variables in several MR studies [5, 6, 30] and including the PREVEND cohort [4] which was employed for these analyses. In addition, the power to detect an effect similar to the observational estimates was low and larger samples (with higher number of cases) are required. In addition, the availability of large samples for genome-wide association study, such as UK Biobank [51], will increase the discovery of genetic variants, which will in turn provide better instruments for MR analyses.

Conclusion

Elevated circulating total bilirubin was independently associated with decreased risk of new onset NAFLD in a cohort of apparently healthy predominantly Caucasian men and women without pre-existing CVD or NAFLD. Based on MR analysis, there was little evidence to suggest a causal association. The observational association may be driven by biases such as unmeasured confounding and/or reverse causation. However, due to low statistical power, larger-scale investigations are necessary to draw definitive conclusions.