Does Bilirubin Have a Causal Relationship With Vitiligo? A Mendelian Randomization Study and Bioinformatics Analysis

Introduction

Vitiligo affects 0.1% to 2% of the global population, presenting as an acquired pigmentary disorder marked by the loss of epidermal melanocytes.^1,2 Its development is associated with genetic predisposition, autoimmune processes, oxidative stress, and the production of inflammatory cytokines.³ Oxidative stress is regarded as a key factor in vitiligo onset.⁴ The accumulation of reactive oxygen species (ROS) disrupts melanocyte integrity, compromising their structural and functional stability, and damaging DNA, lipids, and proteins.⁵ This disturbance triggers various cell death mechanisms, including apoptosis, autophagy, autophagic cell death, and ferroptosis.⁶ Vitiligo manifests in two primary forms: segmental vitiligo (SV) and non-segmental vitiligo (NSV). SV is characterized by localized depigmented patches confined to specific areas, such as the arms or legs, while NSV presents as asymmetrical depigmented lesions distributed across multiple body regions.⁷ Generalized vitiligo (GV), a subtype of NSV, typically involves multifocal and bilateral lesions and is more strongly associated with genetic predisposition than other subtypes.^8–10 Despite considerable progress in understanding its pathogenesis, the exact molecular mechanisms of vitiligo remain unknown, warranting further research.

Bilirubin, the end product of heme metabolism in mammals, is traditionally considered a lipophilic waste that requires elimination from the body due to its potential cytotoxicity.¹¹ While high concentrations can induce cytotoxic effects, including brain damage, recent studies in both in vivo and in vitro models have demonstrated that bilirubin, at physiological levels, possesses significant antioxidant properties.^12,13 It effectively scavenges ROS and inhibits nicotinamide adenine dinucleotide phosphate (NADPH) oxidase activity, thereby reducing oxidative stress and protecting against damage induced by pro-oxidants.^14,15 In addition to its role as an antioxidant in conditions like diabetes, obesity, and cardiovascular diseases, bilirubin levels have been found to decrease in certain skin disorders characterized by oxidative stress. Studies have reported lower bilirubin concentrations in conditions such as pemphigus vulgaris, psoriasis vulgaris, and rosacea.^16,17 Notably, D. Türkmen et al observed significantly reduced serum levels of total bilirubin (TBIL), direct bilirubin (DBIL), and indirect bilirubin (IBIL) in vitiligo patients compared to healthy controls.¹⁸ Notably, Zhang et al observed a marked reduction in serum TBIL levels among vitiligo patients, likely linked to reduced heme oxygenase-1 (HO-1) activity in regulatory T cells (Tregs).¹⁹ This observation highlights bilirubin’s potential involvement in vitiligo pathogenesis. Additionally, an AI-driven fusion of multisource data has identified bilirubin as a significant diagnostic marker for vitiligo.²⁰ However, the exact relationship between bilirubin and vitiligo remains unclear, particularly regarding its role in the pathogenesis of vitiligo via oxidative stress mechanisms. Thus, further research is needed to elucidate bilirubin’s function in vitiligo, enhancing the understanding of its pathological processes and informing potential therapeutic approaches.

To bridge the existing knowledge gap, a Mendelian randomization (MR) analysis was performed to examine the causal relationship between exposure and outcome. This approach leverages genetic variation as an instrumental variable (IV), representing a novel methodology in epidemiological research.^21,22 Using genome-wide association study (GWAS) data, the study aimed to clarify the causal connection between bilirubin and vitiligo. Additionally, a comprehensive bioinformatics analysis, incorporating Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses, was conducted to investigate the molecular mechanisms involved. The findings offer new insights into the bidirectional causality between bilirubin and vitiligo, while also providing a detailed exploration of the underlying mechanisms. This work substantially enhances diagnostic strategies and therapeutic approaches for vitiligo.

Materials and Methods

Study Design

Figure 1 provided a schematic overview of the study design.

Figure 1 Study flow diagram: The research encompassed two primary components: (A) a bidirectional two-sample Mendelian randomization analysis aimed at investigating the relationship between three types of bilirubin (total, direct, and indirect) and vitiligo; (B) a bioinformatics analysis investigating the association between bilirubin and three subtypes of vitiligo.

Data Sources and Summary

The GWAS summary data for vitiligo and three types of bilirubin were retrieved from the Integrative Epidemiology Unit (IEU) OpenGWAS database (https://gwas.mrcieu.ac.uk/). The vitiligo dataset (finn-b-L12_VITILIGO) included 131 cases and 207,482 controls, with 16,380,442 single nucleotide polymorphisms (SNPs) analyzed. The TBIL dataset (ukb-d-30840_irnt) and the DBIL dataset (ukb-d-30660_irnt) comprised 13,585,986 and 13,584,679 SNPs, respectively. The IBIL dataset (ukb-e-recode1_CSA) consisted of 6972 samples and 9,810,003 SNPs.

Data Pre-Processing

The selection of appropriate IVs is fundamental in MR analysis.²³ This method employs genetic variants associated with the exposure to assess causal relationships between the exposure and the outcome. Identifying suitable IVs requires adherence to three core assumptions: (1) a strong association must exist between the IVs and the exposure; (2) IVs must remain unaffected by confounding variables; and (3) the effect of the IVs on the outcome should be mediated solely through their influence on the exposure, with no alternative pathways. Exposure factors were extracted and IVs filtered using the “extract_instruments” function in the R package “TwoSampleMR” (P < 5×10⁻⁸),²⁴ followed by the selection of IVs significantly associated with exposure factors. To minimize duplication and avoid biased estimates of causal effects, independent SNPs were clumped (clumP = TRUE) at a linkage disequilibrium (LD) threshold of r² = 0.001 within a genomic window of 10,000 kb.²⁵ SNPs with a minimum minor allele frequency (MAF) > 0.01 were retained to ensure data quality. The effect estimates were subsequently standardized using the “harmonise_data” function from the “TwoSampleMR” package to ensure consistency in effect sizes. In the forward MR analysis, the exposure was the three bilirubin types, with vitiligo as the outcome. In the reverse MR analysis, these roles were reversed, with vitiligo as the exposure and bilirubin as the outcome.

Bidirectional MR Analysis

The bidirectional MR analysis was conducted using five algorithms in conjunction with the “mr” function: MR Egger, Weighted Median, Inverse Variance Weighted (IVW), Simple Mode, and Weighted Mode. The IVW method served as the primary analytical approach for deriving the MR results. Causal relationships with significant statistical evidence (P < 0.05) were visualized using scatter plots, forest plots, and funnel plots. Cochran’s Q statistic was employed to assess heterogeneity among individual SNPs in causal effect distributions. Sensitivity analyses were performed through three complementary methods: (1) Heterogeneity assessment, where a Q value > 0.05 indicated no heterogeneity; (2) Evaluation of horizontal pleiotropy, where a p-value > 0.05 suggested absence of pleiotropy; and (3) Leave-one-out (LOO) analysis, where the meta-effect of the remaining SNPs was computed after sequential exclusion of each SNP, confirming the reliability of MR results without the influence of SNPs showing substantial bias. This methodological approach ensures analytical robustness and accounts for potential biases in causal inference. Supplementary Methods contain the full implementation details and code used in this analysis.

Functional Annotation

The target SNPs were analyzed using the SNP database from the eQTLGen Consortium (https://eqtlgen.orgphase1.html), which facilitated the identification of genes linked to these SNPs. By inputting the SNP names, corresponding cis-eQTLs were retrieved (cis-regulated, https://eqtlgen.org/cis-eqtls.html). Cis-eQTLs refer to genetic variants situated near the genes they influence, typically within a 1 Mb window upstream or downstream of the gene.

Data Sources for Bioinformatics Analysis

Bilirubin-related genes were extracted from the GeneCards database (version 5.10, available at https://www.genecards.org) using the keyword “bilirubin.” Selection was based on a relevance score of 5 or higher, in accordance with established protocols. Two vitiligo blood datasets were screened: GSE80009 (GPL16951), which included four cases of generalized vitiligo, four cases of segmental vitiligo, and four controls; and GSE90880 (GPL8300), which comprised eight non-segmental vitiligo samples and six control samples for single nucleus RNA sequencing.

Bioinformatic Analysis

Differentially expressed genes (DEGs) were identified, followed by GO and KEGG enrichment analyses performed using SangerBox (http://sangerbox.com). A hypergeometric test assessed the over-representation of GO terms and KEGG pathways. Enrichment levels for both GO and KEGG pathways were evaluated through the same test. P-values were adjusted using the Benjamini-Hochberg (BH) method to control the false discovery rate (FDR). A significance threshold of an adjusted p-value below 0.05 was applied to determine the presence of significant enrichment in pathways or categories.

Ethics Statement

This study performed a secondary analysis using publicly accessible aggregate summary data from prior research and public databases. No original data were collected, nor was there any direct interaction with study participants. In accordance with Article 32 of the Measures for Ethical Review of Life Science and Medical Research Involving Human Subjects, issued by China on February 18, 2023, research utilizing legally acquired public or observational data that does not interfere with public behavior, along with studies using anonymized data, was exempt from institutional ethical approval. Ethical approvals for the primary studies were provided in the respective original publications.

Result

Total and DBIL as Protective Factors Were Causally Correlated With Vitiligo

In the forward MR analysis, 73 SNPs were identified as IVs for DBIL, 3 SNPs for IBIL, and 126 SNPs for TBIL, all of which showed no association with vitiligo (Supplementary Table S1). As presented in Table 1, TBIL (P = 0.038) and DBIL (P = 0.013) were causally linked to vitiligo, while no such association was observed for IBIL (P = 0.157). Additionally, TBIL (b = −0.322) and DBIL (b = −0.460) acted as protective factors for vitiligo, as indicated by the negative slopes in the scatter plots (Figure 2A). Forest plot analysis further assessed the diagnostic performance of each SNP associated with these factors, revealing statistically significant combined effects in both the MR Egger and IVW models (Figure 2B). The funnel plot, used to assess the randomness of the analysis, indicated consistency with Mendel’s second law of random assortment (Figure 2C). In conclusion, total and DBIL demonstrated causal relationships with vitiligo.

Table 1 Mendelian Randomization Analysis of Three Types of Bilirubin (Exposure) and Vitiligo (Outcome). Bold Font Highlights Key Results in This Table, Which Include: 1. The Primary Analytical Method (Inverse Variance Weighted, IVW), Essential to Mendelian Randomization Analysis; 2. p-values

Figure 2 (A) scatter plots for the causality between direct bilirubin/total bilirubin and vitiligo. (B) forest plots for the causality between direct bilirubin/total bilirubin and vitiligo. (C) funnel plots for the causality between direct bilirubin/total bilirubin and vitiligo.

Reliability of the Forward MR results Was Demonstrated Through Sensitivity Analysis

Following the MR analysis, a sensitivity analysis was conducted to evaluate the robustness of the results. The IVW method revealed that the Q_pval from the heterogeneity test and the P values from the horizontal pleiotropy test exceeded 0.05 for both TBIL (Q_pval = 0.470, P = 0.080) and DBIL (Q_pval = 0.775, P = 0.831), indicating no evidence of heterogeneity or horizontal pleiotropy (Tables 2 and 3). Additionally, the LOO method, excluding one SNP at a time, showed no significant impact on the model’s outcomes, suggesting that no single SNP exerted a disproportionate influence (Figure 3). In summary, TBIL and DBIL were confirmed as protective factors against the development of vitiligo, with their effects demonstrating reliability.

Table 2 The Q_pval From the Heterogeneity Test for Total Bilirubin and Direct Bilirubin

Table 3 The P values for the Horizontal Pleiotropy Test Concerning Total Bilirubin and Direct Bilirubin

Figure 3 Leave-one-out plots for the causality between direct bilirubin/total bilirubin and vitiligo.

Reverse MR Analysis Showed No Causal Effect of Vitiligo on Direct Bilirubin

Based on the forward MR results, DBIL, a component of TBIL, was identified as the primary factor in the causal onset of vitiligo. A subsequent reverse MR analysis was performed to examine the relationship between vitiligo and DBIL. Following screening, four SNPs not associated with DBIL were selected as IVs for the reverse MR analysis (Supplementary Table S2). Table 4 presented the finding that vitiligo did not influence DBIL concentration alterations, as determined by five distinct methods (P > 0.05). No causal effect of vitiligo on DBIL levels was observed, as indicated by the reverse MR analysis across all methods (P > 0.05).

Table 4 The Reverse Mendelian Randomization Analysis of Vitiligo and Direct Bilirubin

Functional Enrichment Analyses

To investigate the molecular mechanisms underlying the connection between bilirubin and vitiligo pathogenesis, a comprehensive bioinformatics analysis was performed, incorporating differential gene expression and functional enrichment strategies. Gene expression data from the GEO database were used to examine three vitiligo subtypes—segmental vitiligo (SV), non-segmental vitiligo (NSV), and generalized vitiligo (GV). DEGs were identified for each subtype and cross-referenced with bilirubin-related genes from the Genecards database (association score ≥5) to identify common genes that may mediate the relationship between bilirubin and vitiligo. These co-DEGs were subsequently analyzed through GO and KEGG pathway enrichment to explore their biological functions and associated signaling pathways.

In GV (Figure 4A and B), a total of 136 co-DEGs were identified through intersection. Biological Process (BP) analysis revealed enrichment in response to oxidative stress, regulation of cytokine production, leukocyte proliferation, and other related functions. These genes are localized to membrane rafts and the endoplasmic reticulum lumen in Cellular Component (CC). Molecular Function (MF) assessment highlights antioxidant and ATPase-coupled transmembrane transporter activities. KEGG analysis showed significant enrichment of DEGs in the PI3K-Akt, HIF-1 signaling pathways, as well as in Th17, Th1, and Th2 cell differentiation. In SV (Figure 4C and D), 32 co-DEGs were identified through intersection. BP analysis indicated involvement in the acute inflammatory response, response to hypoxia, and the extrinsic apoptotic signaling pathway, among other functions. These DEGs are localized to membrane rafts and peroxisomes in CC. MF-related DEGs include cytokine and interleukin-1 receptor binding. KEGG analysis revealed prominent enrichment in the NF-kappa B pathway, IL-17 signaling, and Th17 cell differentiation. In NSV (Figure 4E and F), 9 co-DEGs were identified through intersection. BP analysis emphasized the regulation of angiogenesis, JAK-STAT receptor signaling, and response to reactive oxygen species. These genes are localized in the cytoplasmic vesicle lumen, secretory granule lumen, and proteasome core complex in CC. MF analysis showed enrichment in cytochrome-c oxidase activity and DNA binding, bending. KEGG analysis identified significant enrichment in the JAK-STAT signaling pathway and steroid hormone biosynthesis. Based on the identified corresponding cis-eQTLs, further analysis was conducted in conjunction with GV, SV, and NSV. The results revealed that GV intersected with 9 DEGs, SV intersected with 1 DEG, while NSV showed no intersection with any DEGs (Figure 4G).

Figure 4 The Venn diagram shows the overlapping DEGs from the intersection of GV, SV, and NSV with healthy individuals and bilirubin-related genes. GO analysis reveals significant biological processes and molecular functions associated with these gene sets, while KEGG pathway analysis highlights key metabolic pathways involved. (A) GV Venn diagram. (B) GV: GO and KEGG enrichment analysis. (C) SV Venn diagram. (D) SV: GO and KEGG enrichment analysis. (E) NSV Venn diagram. (F) NSV: GO and KEGG enrichment analysis. (G) Overlapping genes between eQTLGen SNPs and vitiligo subtypes.

In summary, these results demonstrate distinct but overlapping molecular pathways linking bilirubin to vitiligo pathogenesis across its subtypes. The identified co-DEGs and enriched pathways highlight the role of oxidative stress, immune regulation, and melanocyte survival in vitiligo, offering mechanistic insights into the protective effects of bilirubin and suggesting potential therapeutic targets for further exploration.

Discussion

This study employed MR to investigate the causal relationship between bilirubin and vitiligo. The results demonstrated a significant inverse association between elevated TBIL and DBIL levels and reduced vitiligo risk, suggesting that both TBIL and DBIL act as protective factors against vitiligo. In contrast, reverse MR analysis found no causal effect of vitiligo on bilirubin, supporting a unidirectional relationship. This unidirectional causality emphasizes the role of redox balance in vitiligo pathogenesis, likely mediated by the antioxidant properties of bilirubin.¹¹

Vitiligo pathogenesis results from a complex interaction of genetic predisposition, autoimmune responses, oxidative stress, intrinsic melanocyte defects, and microenvironmental imbalances. Melanocytes are particularly vulnerable to oxidative damage, with excess ROS triggering apoptosis and immune-mediated destruction of these cells.⁶ Recent studies also highlight a strong association between vitiligo and autoimmune diseases, especially autoimmune thyroid disorders, further supporting the idea that vitiligo is a systemic disease rather than merely a localized skin condition.²⁶ Bilirubin has emerged as a potential protective factor due to its antioxidant properties. It scavenges lipid peroxides and reactive oxygen species (ROS), alleviating oxidative stress, while also inhibiting the release of pro-inflammatory cytokines such as TNF-α and IL-6, thereby mitigating immune-mediated damage.²⁷ Although primarily produced in the liver and spleen, bilirubin is transported via the bloodstream to the skin, where it may exert localized effects on melanocytes.^28,29 Recent research also shows that human keratinocytes can synthesize bilirubin in vitro and in 3D human skin models,³⁰ further highlighting its potential role in the skin microenvironment. Bilirubin’s antioxidant properties suggest it may exert localized effects on skin cells, potentially protecting melanocytes from apoptosis in vitiligo patients. This aligns with our MR study, which provides robust evidence for a causal relationship between bilirubin and vitiligo, confirming its protective role through antioxidant properties that mitigate oxidative stress and immune-mediated damage while highlighting its potential systemic and localized therapeutic applications.

Bioinformatics analysis of microarray data from patients with GV and a bilirubin-related gene set revealed 136 DEGs in GV, 32 in SV, and 9 in NSV. Previous studies suggest that all vitiligo subtypes are associated with inflammatory or immune-related mechanisms; however, clinical features and molecular mechanisms vary across subtypes.⁸ GV is marked by a bilateral symmetrical distribution of lesions and exhibits a stronger correlation with immune mechanisms and genetic susceptibility during melanocyte dysfunction.³¹ Consistent with these distinctions, DEGs in GV are enriched in processes such as autophagy, apoptosis, melanocyte biology, protein degradation, and tyrosine metabolism.⁹ Our findings further corroborate these observations. The higher number of co-DEGs identified in GV compared to SV and NSV reflects more pronounced genetic involvement in GV. These results are consistent with our eQTL analysis, highlighting that GV exhibits a stronger genetic susceptibility compared to other subtypes. Furthermore, generalized vitiligo is widely recognized as a systemic autoimmune disease, featuring widespread immune activation and a pronounced relationship with oxidative stress. This heightened oxidative stress is likely intertwined with bilirubin metabolism, given bilirubin’s role as an endogenous antioxidant. These findings support a potential causal relationship between bilirubin and vitiligo, implicating bilirubin metabolism as a key factor in the underlying pathogenic mechanisms of generalized vitiligo.

A critical finding from this study is the role of oxidative stress, a process closely linked to vitiligo pathogenesis, though manifested differently across subtypes. In GV, DEGs are enriched in pathways like “response to oxidative stress” and “antioxidant activity”, reflecting a strong reliance on antioxidant defense mechanisms, including bilirubin metabolism. The enrichment of ATPase-coupled transport activity in GV further supports its role in maintaining intracellular redox balance. By contrast, SV is characterized by pathways such as “response to hypoxia” and “the extrinsic apoptotic signaling pathway”, where hypoxia-driven ROS accumulation disrupts melanocyte homeostasis and triggers apoptosis through external signals. In NSV, enrichment in pathways like “ROS response” suggests a unique oxidative stress profile, potentially driven by reduced cytochrome c oxidase activity, which limits mitochondrial ROS generation but alters ROS signaling.^32,33 These findings support bilirubin’s role as a potent antioxidant, protecting melanocytes by scavenging ROS, inhibiting lipid peroxidation, and maintaining redox balance, highlighting its potential as a key factor in vitiligo pathogenesis and a therapeutic target.

The subtypes of vitiligo and their associated pathological pathways related to bilirubin demonstrate considerable heterogeneity. In GV, DEGs are notably enriched in the PI3K-Akt, HIF-1, and Th cell differentiation pathways. The PI3K-Akt pathway, a key regulator of cell survival and apoptosis, may enable bilirubin to mitigate oxidative damage to melanocytes in vitiligo through modulation of this pathway.³⁴ The enrichment of Th1/Th2/Th17 cell differentiation pathways further suggests that immune dysregulation—specifically Th17-mediated autoimmune responses—plays a significant role in the pathogenesis and in mediating the protective effects of bilirubin in GV. In contrast, the activation of NF-κB and IL-17 signaling pathways in SV likely reflects direct melanocyte damage resulting from a localized inflammatory microenvironment. Additionally, the enrichment of the JAK-STAT pathway in NSV is consistent with recent evidence supporting the clinical effectiveness of JAK inhibitors in vitiligo treatment.³⁵ These observations provide a theoretical framework for future combination therapies involving bilirubin enhancement, such as the use of JAK inhibitors in conjunction with antioxidant strategies, opening new avenues for therapeutic exploration.

This study indicates that moderately elevated serum bilirubin may protect against vitiligo. However, the potential for excessively high bilirubin concentrations to aggravate oxidative stress and exacerbate vitiligo remains unexamined. While the study offers preliminary insights into the causal relationship and pathway characteristics between bilirubin and vitiligo, the relatively small sample size may limit the robustness of subtype differences. Additionally, the precise mechanisms underlying the interaction between DEGs and bilirubin—such as whether bilirubin directly modulates the PI3K-Akt or JAK-STAT pathways—warrant further experimental investigation. Future studies should aim to expand the cohort, incorporate diverse ethnic groups, and evaluate the potential synergistic effects of other antioxidants (eg, glutathione) alongside bilirubin to fully assess their therapeutic potential in diseases related to oxidative stress.

Recent studies highlight the therapeutic potential of topical bilirubin in various skin conditions. Topical bilirubin alleviated inflammation in psoriasis by modulating MMP9, the MAPK pathway, and cytokines.³⁶ Bilirubin nanoparticles promoted wound healing by balancing cytokines, enhancing angiogenesis, and tissue remodeling,³⁷ while a bilirubin-DFO combination showed antioxidant and pro-angiogenic effects in diabetic wounds.³⁸ These findings lay the groundwork for exploring bilirubin’s potential as a therapeutic option for vitiligo with promise for clinical application.

Strengths and Limitations

The application of MR in this study represents a significant strength, effectively addressing reverse causality and minimizing the influence of confounding variables often present in observational studies, thereby offering clearer insights into the causal links between bilirubin and vitiligo. The incorporation of robust GWAS data further enhances the analysis by providing a reliable basis for identifying IVs associated with bilirubin. In addition, comprehensive bioinformatics analyses, including enrichment studies and core gene identification, are performed to elucidate the molecular mechanisms underlying the relationship between bilirubin and vitiligo.

However, this study is subject to several limitations. First, the GWAS data predominantly represent European populations, while the GEO data primarily originate from Asia and the United States, which may limit the applicability of the findings to other ethnic groups. While this study provides valuable insights, certain limitations should be acknowledged. The eQTL analysis was influenced by the current database availability, which may have constrained deeper analyses. Additionally, the integration of SNP data related to associated diseases presented challenges in fully excluding SNPs potentially linked to liver function, which could influence bilirubin synthesis. Future studies with larger, more diverse cohorts, improved eQTL datasets, detailed liver function data, and skin lesion analyses are needed to validate these findings, provide direct insights into bilirubin’s role in vitiligo, and elucidate the underlying pathogenic mechanisms.

Conclusion

This study presents the first genetic evidence supporting bilirubin’s protective role against vitiligo, with significant inverse associations observed for both total and DBIL. This protective effect is likely mediated by bilirubin’s antioxidant properties, which mitigate melanocyte oxidative stress, a key driver of vitiligo pathogenesis. MR analyses confirm the unidirectional nature of this relationship, with no causal influence of vitiligo on bilirubin. Bioinformatics analysis further validates these findings, identifying enrichment in oxidative stress pathways and differential gene expression patterns across vitiligo subtypes, aligning with bilirubin’s proposed mechanism of action.

Although these results suggest the therapeutic potential of targeting bilirubin pathways in conjunction with immunomodulation, limitations in sample diversity and the need for experimental validation of molecular mechanisms must be acknowledged. Future research should focus on functional studies to translate these genetic insights into clinical applications.