Development and external validation of a mitophagy-related diagnostic model for biliary atresia based on cellular infiltration patterns

Introduction

Biliary atresia (BA) is a severe biliary tract disease characterized by obstruction of the extrahepatic bile ducts, leading to early-onset pathological jaundice and liver failure in infants. The incidence of BA is relatively high in Asia and the Pacific region (1). Untimely diagnosis can result in the destruction of both extrahepatic and intrahepatic bile ducts, disrupting bile flow and causing rapid, severe, and irreversible liver fibrosis (2). Immediate surgical intervention is necessary post-diagnosis, as untreated BA can progress rapidly to end-stage liver cirrhosis within 1 year (3). Although the complete etiology of BA remains incompletely understood, some hypotheses suggest potential associations with immune function (4) and mitochondrial function (5). Molecular biomarkers can accurately capture key biological alterations during disease progression, thereby enabling more precise and efficient diagnosis. In recent years, numerous studies have explored the diagnostic potential of serum indicators (6), bile acid levels (7), and immune-related features (8) in BA; however, their overall performance remains limited, highlighting an urgent need for reliable molecular markers to improve early detection.

Numerous studies have highlighted the enhanced efficacy of autoimmune-modifying therapies in treating BA (9). However, a more comprehensive understanding of the underlying mechanisms of BA is crucial not only to substantiate the effectiveness of immunotherapies but also to pave the way for the exploration of novel treatment strategies. Mitochondria play a pivotal role as key contributors to the cascade of innate immune signaling (10). Mitophagy, a process to remove impaired mitochondria, has been shown to alleviate oxidative damage (11), correlating with liver inflammation and fibrosis levels (12). These observations potentially hint at underlying therapeutic mechanisms for BA. Currently, there is a limited in-depth exploration of the molecular mechanisms involving mitophagy-related genes in BA. Thus, the impact of mitophagy-related gene expression on the autoimmune state (13) may offer crucial insights into the fundamental pathogenic mechanisms and potential drug targeting strategies for BA.

This study aimed to develop and externally validate a multivariable diagnostic prediction model for BA based on mitophagy-related genes and cellular infiltration signatures. In this research, we meticulously analyzed the liver transcriptome of 31 infants with BA, comparing it with data from 20 control liver tissues (Ctrl) with hepatoblastoma. The analysis process is shown in Figure 1. Based on various analyses, we propose the Biliary Atresia Cell Infiltration Mitophagy Diagnostic Gene Model (BA-CIMDGM). Furthermore, drug delivery strategies for BA-CIMDGM were determined based on simulations using AlphaFold3 and Autodock. We present this article in accordance with the TRIPOD reporting checklist (available at https://tp.amegroups.com/article/view/10.21037/tp-2025-1-864/rc).

Figure 1 Research concept and design. BA, biliary atresia; DEGs, differential expression genes; SVM, support vector machine; WGCNA, weighted gene co-expression network analysis.

Methods

Data collection and quality control of BA

The training dataset included 51 participants from Beijing Children’s Hospital, comprising 31 infants with BA and 20 control subjects with hepatoblastoma (Ctrl). The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. The study was approved by the Ethics Committee of Beijing Children’s Hospital Affiliated with Capital Medical University (No. 2023-E-126-Y) and informed consent was obtained from all the parents of each patient. The outcome was defined as a confirmed clinical diagnosis of BA established by pediatric surgeons according to standard diagnostic procedures, including cholangiography, intraoperative evaluation, and histopathological examination. Because this was a retrospective analysis of previously collected clinical data, the assessment of the outcome was not blinded to predictor information.

Liver samples underwent transcriptome sequencing using Illumina RNA-seq (14). Additionally, validation involved collaborative exploration and validation with Gene Expression Omnibus (GEO) datasets related to BA, including GSE122340, GSE15235, and GSE46995 collections, consisting of 296 BA patients and 24 age-matched Ctrl subjects. Expression data were log2 normalized, and averaging yielded final gene expression values. Single-cell (15) (GSE176189) sequencing data for BA field were integrated to investigate cell-level gene expression features.

Physiological phenotypes included aspartate aminotransferase (AST), alanine aminotransferase (ALT), γ-glutamyl transpeptidase (GGT), total bilirubin (TBIL), direct bilirubin (DBIL), indirect bilirubin (IBIL), total bile acid (TBA), liver cirrhosis, and jaundice treatment. To identify mitophagy-related genes (MRGs) involved in BA, we conducted a comprehensive search on the GeneCards database (https://www.genecards.org/) on January 3, 2024. Using the keyword ‘mitophagy’, we initially identified a total of 4,910 potential genes. To ensure the functional significance and biological relevance of the candidates, we applied a Relevance Score threshold of >1. This filtering process resulted in a refined set of 1,672 mitophagy-related genes, which were used for subsequent analyses. The complete list of these genes and their respective relevance scores are provided in table available at https://cdn.amegroups.cn/static/public/tp-2025-1-864-1.xlsx.

Gene expression characteristics and cell infiltration in patients with BA

Transcriptomic landscape analysis of BA and Ctrl

Utilizing a two-sided t-test based on the gene expression data from 31 BA and 20 Ctrl samples, differential expression genes (DEGs) were determined (P value <0.05, |logfold change (FC)| >0.585). Additionally, the Bray-Curtis distance was computed using the Vegdist function to assess the overall transcriptomic features between each pair of samples. Subsequently, principal coordinates analysis (PCoA) clustering was performed to explore the variation in phenotypic features among samples, with continuous physiological phenotypes serving as classification criteria. All these analyses were conducted in the R 4.3.2 programming language, employing packages such as limma 3.9.19 (16) and ggpubr (v 0.4.0).

Differences in cellular infiltration between BA and Ctrl

Cellular infiltration status, based on MCPcounter, EPIC, xCell, CIBERSORT, IPS, ESTIMATE, TIMER, and quanTIseq methods, calculated the distribution proportions of 130 cell types. The differences in cell proportions between BA and Ctrl populations were compared using the Mann-Whitney U test (P value <0.05). Additionally, Spearman correlation analysis was performed to explore the correlation between single-gene expression and overall cellular infiltration across all samples (BA and Ctrl) (P<0.05, |Coef.| >0.3). All these analyses were conducted in the R 4.3.2 language.

Identification of co-expression network of BA

Co-expressed network of mitophagy genes

A weighted gene co-expression network analysis (WGCNA) was conducted based on the expression data of 1,619 mitophagy genes (available online: https://cdn.amegroups.cn/static/public/tp-2025-1-864-1.xlsx, excluding 53 genes with expression values exceeding 80% being 0) (17). WGCNA is a method used to identify sets of genes that exhibit highly coordinated changes at the transcriptional level. In this study, six distinct effective clusters were identified. During WGCNA, genes within the top 75% of median absolute deviation (MAD >0.01, maxBlockSize =5,000, R² =0.85) were excluded to determine the soft threshold, and outlier samples were also removed. All these analyses were conducted in R 4.3.2 language using the WGCNA package (v 1.71).

Co-expressed network related to cellular infiltration

A Spearman correlation test was performed based on the differences in cellular infiltration between BA and Ctrl, along with the six effective clusters from the WGCNA network. This was done to identify characteristic networks associated with BA (P value <0.05, |Coef.| >0.2). Finally, the feature genes in the network were filtered based on transcriptional differences between BA and Ctrl (P value <0.05, |logFC| >0.585).

Feature gene selection and model construction of mitophagy in BA

BA immune-related feature gene selection

The study utilized least absolute shrinkage and selection operator (LASSO) and support vector machine (SVM) algorithms to comprehensively screen transcriptional data for significant biomarkers, leveraging model coefficient sparsity and optimal hyperplane selection, respectively. The intersection of their results identified the optimal distinguishing features between BA and Ctrl. The LASSO regression utilized a 20-fold cross-validation (binomial) for binary feature selection. Additionally, SVM employed the svmRadial method for screening. The Venn diagram method was applied to extract common feature genes identified by both LASSO and SVM algorithms. Multivariate regression modeling was employed to establish the immunodiagnostic model for BA using feature genes, excluding genes with a P value less than 0.5. The validation process (GEO) utilized the model, along with gene expression values estimated by the model, to calculate the area under the curve (AUC), thereby investigating the diagnostic capabilities of both multi-gene and single-gene approaches. The above analysis was performed in R (v 4.3.2) language.

External validation of the BA immune diagnostic model genes

Due to significant differences in sample sizes between BA and Ctrl across the three GEO datasets, we employed a merged approach for transcriptome differential validation. The ComBat method (18) was employed to eliminate transcription batch effects in three GEO datasets (GSE122340, GSE15235, and GSE46995) (Figure S1). Furthermore, a two-sided T-test was employed to explore the consistency of gene expression differences between the validation and training datasets (P value <0.05). This analysis aimed to assess the consistency of gene expression across different datasets and ascertain the robustness of the BA immune diagnostic model. The above analysis was performed in R (v 4.3.2) language.

Single-cell transcriptome (scRNA-seq) model of gene immune cell expression states in BA

scRNA-seq data from nine BA patients and three control tissues were obtained from Luo et al. (15). A reanalysis of the data was conducted using the Seurat software package (19). The analysis aimed to comprehensively characterize the expression profiles of feature genes at the single-cell level. Seurat object inclusion criteria were to retain genes with assay counts between 200 and 5,000 and unique molecule (UMI) counts below 50,000. To address batch effects, the ‘FindIntegrationAnchors’ function was applied. Subsequently, the ‘FindNeighbors’ and ‘FindClusters’ functions (dim =10, resolution =0.5) were used to cluster individual cells into distinct subgroups. Marker genes for each cell subgroup were identified using the ‘FindAllMarkers’ function, associating a cluster with logFC >1. The ACT database (xteam.xbio.top/ACT/) was utilized for annotating cell types in liver tissues. Following annotation, the ‘FindMarkers’ function was employed to calculate significant differences between BA and Ctrl across various cell subgroups (P value adj <0.05). This analysis aimed to determine if there were immune cell expression differences between BA and Ctrl in the model genes. The above analysis was performed in R 4.3.2 language.

Immunohistochemistry staining assay

Liver samples taken from BA patients or adjacent non-tumor liver tissues taken from hepatoblastoma patients who underwent surgery and were fixed in 4% paraformaldehyde and embedded in Paraffin compound. Paraffin-embedded liver sections were firstly dewaxed and dehydrated, subsequently underwent antigen retrieval and inactive endogenous peroxidase procedures, and then were blocked in 5% BSA. Next, sections were incubated overnight at 4 ℃ with the following primary antibodies: anti-human ALDH2 (HUABIO, Hangzhou, China; M1509-1, 1:800), NEDD4L (Proteintech, Wuhan, China; 13690-1-AP, 1:800), and IGF2BP3 (Proteintech, Wuhan, China; 14642-1-AP, 1:800). Then, sections were incubated with horseradish peroxidase (HRP)-labeled secondary antibody (Abcam, Cambridge, UK; ab6721) for 30 min at room temperature followed by detection using DAB substrate and hematoxylin.

Drug screening

Firstly, we screened drugs associated with the three target genes based on the CTD database and extracted their intersection. Then, based on ‘Interaction Actions’, we assessed the drugs’ impact on gene expression (mRNA) and excluded uncertain effects (removing entries labeled ‘affects’), thus identifying the drugs meeting the criteria.

Based on the drug’s structural properties, we explored their potential competitive inhibition effects on the protein

IGF2BP3 and NEDD4L’s chemical structures were simulated using AlphaFold3 (20). The structures of small molecule drugs were sourced from the PubChem database. Subsequently, Autodock (grid and docking) semi-flexible simulations were conducted to determine the active interaction characteristics of the target drugs (with proteins considered rigid and small molecule drugs flexible). The grid box (126 Å × 126 Å × 126 Å, spacing =1) was centered at (IGF2BP3: −1.7758, 9.343, 0.997; NEDD4L: −2.517, 8.64, −1.98) Å (globally for Osteopontin protein). The Docking Genetic Algorithm was configured with 100 Number of GA Runs, Maximum Number of generations set to 27,000, and Maximum Number of evals set to 2,500,000. Conformations were output based on energy (affinity energy/(kcal/mol) <−2 kcal/mol) and interaction atoms. Docking results were visualized using PyMol.

Statistical analysis

All statistical analyses were performed using R software (version 4.3.2). Differential gene expression between BA and Ctrl groups was assessed using two-sided Student’s t-tests (P<0.05, |log₂FC| >0.585). Global transcriptomic differences were evaluated using Bray-Curtis distance and PCoA. Differences in immune cell infiltration were compared using the Mann-Whitney U test, and associations between gene expression and immune infiltration were assessed using Spearman correlation analysis (P<0.05). WGCNA was used to identify mitophagy-related gene modules associated with BA. Feature selection was performed using LASSO regression with 20-fold cross-validation and SVM algorithms. Multivariate logistic regression models were constructed, and diagnostic performance was evaluated using receiver operating characteristic (ROC) curve analysis and AUC. Batch effects in external GEO datasets were corrected using the ComBat method. Single-cell RNA sequencing differential expression was assessed using the Wilcoxon rank-sum test. All tests were two-sided, with P<0.05 considered statistically significant.

Results

Differences in transcriptomic landscape and cell infiltration

This study recruited a total of 51 volunteers (BA =31, Ctrl =20). Based on the differential comparison of transcriptomes, a total of 736 significantly upregulated and 1,472 significantly downregulated DEGs were identified (P<0.05, |logFC| >0.585), respectively. The PCoA of the transcriptome revealed a 65.8% difference in transcriptional features between BA and Ctrl (Figure 2A, PCo1 explained variance), highlighting diverse phenotypic states between BA and Ctrl (Figure 2B).

Figure 2 Differences in transcriptional characteristics between BA and Ctrl. (A) Discrepancies in the immune landscape within the BA population. (B) Phenotypic differences. ***, P<0.001. ALT, alanine aminotransferase; AST, aspartate aminotransferase; BA, biliary atresia; Ctrl, control; DBIL, direct bilirubin; IBIL, indirect bilirubin; TBA, total bile acid; TBIL, total bilirubin; WGCNA, weighted gene co-expression network analysis.

Cell infiltration states showed a significant association with BA (21) and based on the infiltration status of 130 cell types, we identified 13 cell types exhibiting significant differences between BA and Ctrl (Table S1, P<0.05/130). Notably, Tregs, macrophages, mast cells, and other cells displayed elevated infiltration levels in BA (Figure 3A).

Figure 3 WGCNA and cell infiltration status select co-expressed mitophagy gene set. (A) Analysis of cell infiltration. (B) (a) Evaluation of scale independence in WGCNA. (b) Assessment of mean connectivity in WGCNA. (c) Cluster dendrogram in WGCNA. (d) Heatmap displaying eigengene adjacency. (e) Relationships between module and cell infiltration in WGCNA. (C) Twenty-nine genes co-expressed in BA-infiltrated WGCNA clusters and were significantly different between BA and Ctrl. BA, biliary atresia; WGCNA, weighted gene co-expression network analysis.

Identification of co-expression network of BA

To explore the co-expression status of mitophagy genes in BA patients and establish a co-expression network related to cell infiltration states, we identified six clustering modules from 1,619 genes in 45 samples (Table S2, R²=0.85, Power =18, Figure 3B, a-c). Clustering based on module scores revealed 884 genes with similar co-expression patterns in Cluster Red, turquoise, and yellow (Figure 3B, d). Significant associations were found between 13 cell types and three modules (Figure 3B, e; P<0.01, |Coeff| >0.2). Subsequently, we identified 29 mitophagy genes with transcriptional differences in three cell infiltration-related modules between BA and Ctrl (Figure 3C, Table 1).

Establishment and validation of the BA-CIMDGM

Using LASSO and SVM on 29 differentially expressed, cell infiltration-related genes, we identified 10 and 11 genes associated with BA diagnosis (Table S3, Figure 4A, a,b), respectively. This intersection prominently highlights six cell infiltration-related genes significantly differing between BA and Ctrl (Figure 4A, c). Using a multivariate regression model, we re-evaluated the selectivity of these six genes for BA and Ctrl, with four genes demonstrating significant feature levels in the overall model, which were employed to establish the BA-CIMDGM (P<0.05, Figure 4A, d). Testing the diagnostic capability of the four genes for BA in two validation datasets with BA and Ctrl, ROC analysis demonstrated strong predictive diagnostic levels: GSE46995 [AUC =1.00; 95% confidence interval (CI): 1.00, 1.00] and GSE122340 (AUC =0.826; 95% CI: 0.749, 0.923) (Figure 4A, e).

Figure 4 WGCNA and cell infiltration status select co-expressed mitophagy gene set. (A) (a) Feature selection using LASSO. (b) Feature selection using SVM. (c) Common features identified by LASSO and SVM. (d) Establishment of the BA-CIMDGM. (e) The BA-CIMDGM model was validated in the GEO dataset, obtaining the AUC. (B) (a) Validation of model gene expression differences in GEO datasets. (b) The expression differences between BA and Ctrl. (C) Correlation between feature genes and cell infiltration. *, P<0.05; ***, P<0.001. AUC, area under the curve; BA, biliary atresia; BA-CIMDGM, Biliary Atresia Cell Infiltration Mitophagy Diagnostic Gene Model; CI, confidence interval; Ctrl, control; GEO, Gene Expression Omnibus; LASSO, least absolute shrinkage and selection operator; RMSE, root mean square error; SE, standard error; SVM, support vector machine; WGCNA, weighted gene co-expression network analysis.

To validate gene expression differences in the GEO dataset, we merged data from three GEO datasets due to limited control sample sizes. Results showed significant differences in three genes (IGF2BP3, NEDD4L, and ALDH2; Figure 4B, a). Their expression direction was consistent with the training dataset (Figure 4B, b). There is no significant difference in PPIB. While all four genes exhibit notable diagnostic features for BA, post-expression differential validation revealed that only three genes were independently replicated, serving as crucial molecular markers for BA diagnosis. Immunohistochemistry results of the liver tissues also supported these findings. Compared with control group, ALDH2 (a mitochondrial enzyme critical for acetaldehyde detoxification) nearly located at hepatocytes and showed attenuated staining intensity in the BA group, indicative of abnormal liver function. IGF2BP3 (an RNA-binding protein that stabilizes target mRNAs and promotes their translation) exhibited more punctate staining patterns in the BA group, suggesting enhanced stress and autophagy. NEDD4L (an E3 ubiquitin ligase that tags specific protein substrates with ubiquitin for proteasomal degradation) displays accumulated in the BA group, intense granular staining, indicating impaired mitophagy process. It is noteworthy that IGF2BP3 and NEDD4L are highly expressed in the infiltrating immune cells of the bile ducts and portal areas in the BA group, suggesting they may regulate the function of immune cells, leading to damage to the bile ducts (Figures 4C,5A).

Figure 5 Single cell sequencing and immunohistochemistry validate key molecular markers for BA. (A) Representative immunohistochemistry images validating the protein expression of mitophagy-related genes (IGF2BP3, NEDD4L, and ALDH2) in liver tissues from BA patients and controls (scale bar =100 µm). Ctrl images are from non-tumor liver tissue adjacent to hepatoblastoma specimens. Positive staining appears as brown discoloration and infiltrating immune cells were indicated by red arrows. (B) Differential analysis of single-cell transcription. (a) Single-cell’s cell classification landscape in all samples. (b) The single-cell expression differences of IGF2BP3, NEDD4L, PPIB, and ALDH2 between BA and Ctrl groups. BA, biliary atresia; Ctrl, control; FC, fold change.

Verification of immune cell expression differences in the scRNA-seq of BA

To characterize the hepatic distribution patterns of IGF2BP3, NEDD4L, PPIB, and ALDH2 in BA, we analyzed single-cell transcriptomic data from nine BA patients and three controls. The results revealed marked expression differences of these four genes across multiple immune cell populations (Figure 5B, Table S4). Specifically, IGF2BP3 was significantly upregulated in CD8⁺ T cells, Kupffer cells, and memory B cells from BA patients, whereas PPIB showed a notable downregulation in Kupffer cells. In contrast, NEDD4L exhibited significantly reduced expression in Kupffer cells, and ALDH2 was markedly downregulated in endothelial cells. Correlation analyses further demonstrated strong associations between the expression of these genes and immune cell infiltration in BA, most prominently with Tregs (P=9.84×10⁻⁸, Coef. =0.728; Table S5, Figure 4C). In addition, all three genes showed significant correlations with macrophage abundance.

Drug screening and prediction

To explore drugs that effectively modulate all three genes in the model, we searched the CTD database and found 146, 223, and 573 drugs associated with IGF2BP3, NEDD4L, and ALDH2, respectively. Among them, 21 drugs may affect these three genes. In differential analysis, we found that elevation of IGF2BP3 and NEDD4L may lead to BA, while downregulation of ALDH2 may be a protective factor against BA (Table S6). Based on these drugs, we identified Cyclosporine that may decrease the expression of IGF2BP3 and NEDD4L while increasing ALDH2 mRNA expression. Subsequently, at the protein level, we also tested whether these drugs competitively inhibit IGF2BP3 and NEDD4L proteins. We found that Cyclosporine exhibits significant interactions with IGF2BP3 and NEDD4L. Notably, Cyclosporine interacts with the G31 amino acid residue of IGF2BP3 (A =2.8 pm, energy =−10.26 kcal/mol, Figure 6A-6C), potentially inhibiting the RRM 1 structural region of the protein. Additionally, cyclosporine shows even stronger interactions with NEDD4L, forming hydrogen bonds with the E202 and H214 amino acid residues (A =2.7 pm and 2.5 pm, energy =−10.59 kcal/mol, Figure 6D-6F), which may inhibit the WW 1 structural region of the NEDD4L protein. In conclusion, our results suggest that Cyclosporine could serve as potential therapeutic agents for BA.

Figure 6 The interaction between drugs and IGF2BP3 and NEDD4L proteins. (A-C) The interaction between cyclosporine and IGF2BP3 protein. (D-F) The interaction between cyclosporine and NEDD4L protein. Blue represents drug molecules, with red indicating oxygen atoms within the molecules. Green represents amino acid residues of the protein interacting with the drug molecules. The yellow background indicates the average interaction energy between the protein and drugs (average of 10 simulations).

Discussion

The study systematically evaluated BA’s cellular infiltration using 130 cell features. Notably, macrophages, Treg cells, mast cells, eosinophils, and granulocyte-macrophage progenitor cells showed increased infiltration. At the gene level, IGF2BP3, NEDD4L, and ALDH2 were identified as key mitophagy genes linked to BA’s immune status. These genes consistently outperformed others in BA diagnosis. In immune analysis and single-cell clustering, the genes showed significant associations with cells like macrophages and memory B cells. Immunohistochemistry confirmed altered expression of these genes in BA tissues. In summary, our findings indicate a significant association between three genes and the cellular infiltration status in the progression of BA. Moreover, we conducted drug screening targeting the key genes in BA-CIMDGM.

The genes incorporated in the BA-CIMDGM model might cooperatively drive a pathological cascade in BA, extending from immune-mediated injury to progressive fibrosis (22) (Figure 7). IGF2BP3, a key RNA-binding protein and N6-methyladenosine (m6A) reader, is selectively overexpressed in infiltrating CD8⁺ T cells and Kupffer cells in BA livers. By stabilizing pro-inflammatory mRNA transcripts or modulating m6A-dependent pathways, IGF2BP3 may amplify excessive immune activation within the portal microenvironment (23,24). NEDD4L, an E3 ubiquitin ligase that tags specific protein substrates with ubiquitin for proteasomal degradation, is upregulated in macrophage, indicating impaired mitophagy (25). Additionally, downregulation of the mitochondrial detoxifying enzyme ALDH2 in hepatocytes further compromises antioxidant defenses, dysfunctional mitochondria accumulate and generate excessive reactive oxygen species (ROS), promoting biliary epithelial cell (BEC) injury and ductal obliteration, synergistically accelerating hepatic stellate cell activation and the transition from cholestasis to irreversible liver fibrosis (26). Moreover, IGF2BP3 has been implicated in mitophagy regulation via m6A modification and Parkin signaling (27,28), suggesting a potential role in m6A modification mediated mitophagy in immune cells. Collectively, these findings indicate that sustained oxidative stress and feedback-enhanced immune response jointly contribute to BA pathogenesis and may represent potential therapeutic targets.

Figure 7 Potential schematic representation of the BA-CIMDGM genes in the pathogenesis of BA. Upregulated IGF2BP3 in T cells and NEDD4L in macrophages leads to impaired mitophagy and promotes immune activation and cytokines released. And the downregulation of ALDH2 in hepatocytes compromises the mitochondrial antioxidant defense and excessive ROS production, synergistically driving biliary epithelial cell injury and progressive liver fibrosis. BA, biliary atresia; BA-CIMDGM, Biliary Atresia Cell Infiltration Mitophagy Diagnostic Gene Model; ROS, reactive oxygen species.

The analysis of 130 cellular infiltration features suggests potential immune cellular activation in BA. The integration of our model with immune infiltration analyses suggests a plausible pathogenic link between mitophagy and BA progression. Dysregulated mitophagy in key hepatic cells may fuel a vicious cycle of inflammation and fibrosis. For instance, impaired mitophagy in Kupffer cells could lead to the accumulation of damaged, ROS-generating mitochondria, promoting their polarization towards a pro-inflammatory phenotype and driving the secretion of IL-6 and TNF-α (29,30). Concurrently, altered mitophagy in regulatory T cells (Tregs) might compromise their metabolic fitness and suppressive function, contributing to immune dysregulation (31,32). This creates an oxidative stress feedback loop: inflammation increases mitochondrial damage, while defective clearance of these organelles further elevates ROS and pro-inflammatory signals. Ultimately, this microenvironment likely exacerbates the injury and fibro-obliteration of the bile duct epithelium, a hallmark of BA. Notably, there is an observed increase in the infiltration of macrophages M0 and M1 in BA (21). CD8⁺ T cells and Tregs have also been reported to contribute to infiltrative damage in BA (33,34). Additionally, limited studies have indicated the detrimental role of mast cells in BA (33), and reports on eosinophils have shown a significant correlation with BA risk (35). These findings are consistent with our analyzed results in cell infiltration.

Our drug prediction results identified a positive regulatory effect of cyclosporine on the BA-CIMDGM model. Cyclosporine demonstrates not only a significant impact on gene expression but also exhibits distinct interactions with critical active regions of the IGF2BP3 and NEDD4L proteins. These interactions suggest a potential mechanism for the competitive inhibition of protein function. Cyclosporine is an immunosuppressive agent known for its selective and reversible inhibition of T lymphocytes, exhibiting low cytotoxicity (36). Our drug screening nominated cyclosporine A, an inhibitor of cyclophilins, as a candidate therapeutic agent. This finding is mechanistically supported by experimental evidence showing that a unique cyclosporine A derivative, MM284, inhibited inflammatory and fibrotic pathways in a murine model of BA, and its inhibition ameliorates disease progression (37). However, the direct clinical translation of systemic cyclosporine for infant BA is highly challenging. As a potent calcineurin inhibitor, cyclosporine carries significant risks for this vulnerable population, including nephrotoxicity, hepatotoxicity, hypertension, and increased susceptibility to infections and malignancy. While cyclosporine has a defined role in pediatric liver transplantation post-BA (38) and in autoimmune hepatitis (39), its utility as a primary medical therapy for native liver disease in BA remains speculative and is not part of standard care (40,41). Furthermore, its use would require meticulous therapeutic drug monitoring and management of numerous drug interactions. Additionally, cyclosporine has shown promising efficacy in the treatment of conditions such as asthma and primary biliary cirrhosis (42). In summary, our findings suggest a potential therapeutic agent targeting mitophagy for BA.

Furthermore, compared to existing single-molecule biomarkers such as serum MMP-7 or distinct bile acid profiles, our BA-CIMDGM model offers a multi-gene signature that may provide superior diagnostic accuracy by capturing a more comprehensive pathological state. More importantly, while MMP-7 reflects extracellular matrix remodeling and bile acids indicate cholestasis, our model is rooted in mitophagy and its interplay with immune infiltration. This provides deeper biological insight into the underlying cellular pathogenesis of BA, potentially identifying the disease at an earlier stage in the dysfunctional cascade, before extensive fibrosis or cholestasis becomes irreversibly established. Validating this model in prospective cohort could pave the way for a more mechanistic diagnostic tool.

Despite the robust performance and multi-level validation of the BA-CIMDGM model, several limitations warrant consideration. First, concerning the control group, the liver tissues were obtained from the non-tumor periphery of hepatoblastoma specimens. Although histologically normal, the potential systemic or local paracrine effects of the tumor microenvironment on immune infiltration cannot be entirely excluded. Furthermore, the absence of disease controls, such as Alagille syndrome or progressive familial intrahepatic cholestasis, means the model’s specificity in differentiating BA from other cholestatic etiologies remains to be fully established. Second, regarding the model genes, while PPIB was statistically indispensable for the model’s integrative predictive power, it did not exhibit significant differences in certain external cohorts. This inconsistency may stem from platform-specific sensitivities, limited sample sizes, or the inherent transcriptomic heterogeneity across different pathological stages of BA. Finally, from a translational perspective, this study is primarily limited by its retrospective nature and the lack of functional validation in cellular or animal models. Although molecular docking provided preliminary evidence for therapeutic candidates like cyclosporine, these remain bioinformatics-driven hypotheses. Future large-scale, multi-center prospective cohorts and in-depth experimental assays are essential to confirm the model’s clinical utility across diverse ethnicities and to elucidate the underlying molecular mechanisms.

Conclusions

This study clarifies the critical role of mitophagy in BA progression by defining specific liver cell infiltration states. A novel diagnostic model, the BA-CIMDGM, was constructed using four key mitophagy-related genes (NEDD4L, IGF2BP3, ALDH2, and PPIB), demonstrating high diagnostic accuracy (AUC >0.83) in independent datasets. Furthermore, drug screening identified cyclosporine as a promising therapeutic candidate. In conclusion, these findings provide a significant diagnostic tool and a compelling pharmacological strategy, highlighting the BA-CIMDGM model’s potential to advance both the diagnosis and treatment of BA.

Acknowledgments

We are grateful to the patients and the family of the patients who have made this research possible.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tp.amegroups.com/article/view/10.21037/tp-2025-1-864/rc

Data Sharing Statement: Available at https://tp.amegroups.com/article/view/10.21037/tp-2025-1-864/dss

Peer Review File: Available at https://tp.amegroups.com/article/view/10.21037/tp-2025-1-864/prf

Funding: This study received financial support from National Natural Science Foundation of China grants (#82300574, #82400592), Science and Technology Program of the Joint Fund of Scientific Research for the Public Hospitals of Inner Mongolia Academy of Medical Sciences (No. 2024GLLH0202), Beijing Municipal Science & Technology Commission (No. Z211100002921062), Jiangxi (Ganpo) Talent Program-Health Innovation Talent Project (No. gpyc20240208), and Beijing Municipal Natural Science Foundation (#7252043).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tp.amegroups.com/article/view/10.21037/tp-2025-1-864/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. The study was approved by the Ethics Committee of Beijing Children’s Hospital Affiliated with Capital Medical University (No. 2023-E-126-Y) and informed consent was obtained from all the parents of each patient.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Chung PHY, Zheng S, Tam PKH. Biliary atresia: East versus west. Semin Pediatr Surg 2020;29:150950. [Crossref] [PubMed]
2. Zhu JJ, Yang YF, Dong R, et al. Biliatresone: progress in biliary atresia study. World J Pediatr 2023;19:417-24. [Crossref] [PubMed]
3. Tam PKH, Wells RG, Tang CSM, et al. Biliary atresia. Nat Rev Dis Primers 2024;10:47. [Crossref] [PubMed]
4. Ye C, Zhu J, Wang J, et al. Single-cell and spatial transcriptomics reveal the fibrosis-related immune landscape of biliary atresia. Clin Transl Med 2022;12:e1070. [Crossref] [PubMed]
5. Xiao Y, Zhou Y, Lu Y, et al. PHB2 interacts with LC3 and SQSTM1 is required for bile acids-induced mitophagy in cholestatic liver. Cell Death Dis 2018;9:160. [Crossref] [PubMed]
6. Han YJ, Hu SQ, Zhu JH, et al. Accurate prediction of biliary atresia with an integrated model using MMP-7 levels and bile acids. World J Pediatr 2024;20:822-33. [Crossref] [PubMed]
7. Johansson H, Svensson JF, Almström M, et al. Regulation of bile acid metabolism in biliary atresia: reduction of FGF19 by Kasai portoenterostomy and possible relation to early outcome. J Intern Med 2020;287:534-45. [Crossref] [PubMed]
8. Dong R, Jiang J, Zhang S, et al. Development and Validation of Novel Diagnostic Models for Biliary Atresia in a Large Cohort of Chinese Patients. EBioMedicine 2018;34:223-30. [Crossref] [PubMed]
9. Mack CL, Tucker RM, Sokol RJ, et al. Biliary atresia is associated with CD4+ Th1 cell-mediated portal tract inflammation. Pediatr Res 2004;56:79-87. [Crossref] [PubMed]
10. Chen Y, Zhou Z, Min W. Mitochondria, Oxidative Stress and Innate Immunity. Front Physiol 2018;9:1487. [Crossref] [PubMed]
11. Liu L, Zhang W, Liu T, et al. The physiological metabolite α-ketoglutarate ameliorates osteoarthritis by regulating mitophagy and oxidative stress. Redox Biol 2023;62:102663. [Crossref] [PubMed]
12. Parola M, Robino G. Oxidative stress-related molecules and liver fibrosis. J Hepatol 2001;35:297-306. [Crossref] [PubMed]
13. Xu Y, Shen J, Ran Z. Emerging views of mitophagy in immunity and autoimmune diseases. Autophagy 2020;16:3-17. [Crossref] [PubMed]
14. Wang D, Yang S, Zhao Y, et al. Identifying and validating molecular subtypes of biliary atresia using multiple high-throughput data integration analysis. Front Immunol 2022;13:1008246. [Crossref] [PubMed]
15. Xiao MH, Ma D, Wu S, et al. Integrative single-cell and spatial transcriptomic analyses identify a pathogenic cholangiocyte niche and TNFRSF12A as therapeutic target for biliary atresia. Hepatology 2025;81:1146-63. [Crossref] [PubMed]
16. Ritchie ME, Phipson B, Wu D, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 2015;43:e47. [Crossref] [PubMed]
17. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 2008;9:559. [Crossref] [PubMed]
18. Zhang Y, Parmigiani G, Johnson WE. ComBat-seq: batch effect adjustment for RNA-seq count data. NAR Genom Bioinform 2020;2:lqaa078. [Crossref] [PubMed]
19. Butler A, Hoffman P, Smibert P, et al. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat Biotechnol 2018;36:411-20. [Crossref] [PubMed]
20. AlphaFold3 - why did Nature publish it without its code? Nature 2024;629:728. [Crossref] [PubMed]
21. Wang J, Xu Y, Chen Z, et al. Liver Immune Profiling Reveals Pathogenesis and Therapeutics for Biliary Atresia. Cell 2020;183:1867-1883.e26. [Crossref] [PubMed]
22. Hartley JL, Davenport M, Kelly DA. Biliary atresia. Lancet 2009;374:1704-13. [Crossref] [PubMed]
23. Wang J, Du M, Meng L, et al. Integrative analysis implicates the significance of m6A in the liver fibrosis of biliary atresia by regulating THY1. Hepatol Commun 2023;7:e0004. [Crossref] [PubMed]
24. Ma L, Jiang J, Si Q, et al. IGF2BP3 Enhances the Growth of Hepatocellular Carcinoma Tumors by Regulating the Properties of Macrophages and CD8(+) T Cells in the Tumor Microenvironment. J Clin Transl Hepatol 2023;11:1308-20. [Crossref] [PubMed]
25. Jiang B, Bai F, Hu Y, et al. Endothelial major vault protein alleviates vascular remodeling via promoting Parkin-mediated mitophagy. Nat Commun 2025;16:4365. [Crossref] [PubMed]
26. Ma X, Luo Q, Zhu H, et al. Aldehyde dehydrogenase 2 activation ameliorates CCl(4) -induced chronic liver fibrosis in mice by up-regulating Nrf2/HO-1 antioxidant pathway. J Cell Mol Med 2018;22:3965-78. [Crossref] [PubMed]
27. Wang J, Zheng F, Wang D, et al. Regulation of ULK1 by WTAP/IGF2BP3 axis enhances mitophagy and progression in epithelial ovarian cancer. Cell Death Dis 2024;15:97. [Crossref] [PubMed]
28. Sun X, Ye G, Li J, et al. Parkin regulates IGF2BP3 through ubiquitination in the tumourigenesis of cervical cancer. Clin Transl Med 2023;13:e1457. [Crossref] [PubMed]
29. Wu H, Chen G, Wang J, et al. TIM-4 interference in Kupffer cells against CCL4-induced liver fibrosis by mediating Akt1/Mitophagy signalling pathway. Cell Prolif 2020;53:e12731. [Crossref] [PubMed]
30. Xu Y, Tang Y, Lu J, et al. PINK1-mediated mitophagy protects against hepatic ischemia/reperfusion injury by restraining NLRP3 inflammasome activation. Free Radic Biol Med 2020;160:871-86. [Crossref] [PubMed]
31. Li H, Guo Y, Su W, et al. The mitochondria-targeted antioxidant MitoQ ameliorates inorganic arsenic-induced DCs/Th1/Th2/Th17/Treg differentiation partially by activating PINK1-mediated mitophagy in murine liver. Ecotoxicol Environ Saf 2024;277:116350. [Crossref] [PubMed]
32. Xu L, Zhao B, Cheng H, et al. Bergapten enhances mitophagy to regulate intestinal barrier and Th17/Treg balance in mice with Crohn's disease-like colitis via PPARγ/NF-κB signaling pathway. Naunyn Schmiedebergs Arch Pharmacol 2024;397:7589-97. [Crossref] [PubMed]
33. Tucker RM, Feldman AG, Fenner EK, et al. Regulatory T cells inhibit Th1 cell-mediated bile duct injury in murine biliary atresia. J Hepatol 2013;59:790-6. [Crossref] [PubMed]
34. Sun D, Wang D, Jia L, et al. S1P/S1PR4 promotes the differentiation of CD8(+) tissue-resident memory T cells aggravating bile duct injury in biliary atresia. J Hepatol 2026;84:385-98. [Crossref] [PubMed]
35. Yamada Y, Hoshino K, Shimojima N, et al. Idiopathic hypereosinophilic syndrome in a case with ABO-incompatible liver transplantation for biliary atresia complicated by portal vein thrombosis. Pediatr Transplant 2010;14:e49-53. [Crossref] [PubMed]
36. Guada M, Beloqui A, Kumar MN, et al. Reformulating cyclosporine A (CsA): More than just a life cycle management strategy. J Control Release 2016;225:269-82. [Crossref] [PubMed]
37. Iordanskaia T, Malesevic M, Fischer G, et al. Targeting Extracellular Cyclophilins Ameliorates Disease Progression in Experimental Biliary Atresia. Mol Med 2015;21:657-64. [Crossref] [PubMed]
38. Iwatsuki S, Shaw BW Jr, Starzl TE. Liver transplantation for biliary atresia. World J Surg 1984;8:51-6. [Crossref] [PubMed]
39. Franulović OZ, Rajacić N, Lesar T, et al. Cyclosporine induced biochemical remission in childhood autoimmune hepatitis. Coll Antropol 2012;36:973-9.
40. Ando H, Inomata Y, Iwanaka T, et al. Clinical practice guidelines for biliary atresia in Japan: A secondary publication of the abbreviated version translated into English. J Hepatobiliary Pancreat Sci 2021;28:55-61. [Crossref] [PubMed]
41. Calinescu AM, Wildhaber BE. Post-Kasai cholangitis evaluation and management strategies: Review of the literature with insights from the Swiss Biliary Atresia Registry. Semin Pediatr Surg 2024;33:151471. [Crossref] [PubMed]
42. Survase SA, Kagliwal LD, Annapure US, et al. Cyclosporin A--a review on fermentative production, downstream processing and pharmacological applications. Biotechnol Adv 2011;29:418-35. [Crossref] [PubMed]