Abstract Hepatic ischemia-reperfusion injury (HIRI) is a major complication following liver transplantation. Bioinformatic analysis was performed to elucidate the PANoptosis-related molecular mechanisms underlying HIRI. Comprehensive analysis of bulk and single-cell RNA sequencing data from human liver tissue before and after HIRI was performed. Differential expression analysis, weighted gene coexpression analysis, and protein interaction network analysis were used to identify candidate biomarkers. Multiple machine learning methods were utilized to screen for core biomarkers and construct a diagnostic predictive model. Functional and interaction analyses of the genes were also performed. Cellular clustering and annotation, pseudotemporal trajectory, and intercellular communication analyses of HIRI were conducted. Six PANoptosis-associated genes (CEBPB, HSPA1A, HSPA1B, IRF1, SERPINE1, and TNFAIP3) were identified as HIRI-related biomarkers. These biomarkers are regulated by NF-κB and miRNA-155. A nomogram for HIRI prediction based on these biomarkers was constructed and validated. In addition, the heterogeneity and dynamic changes in macrophage subpopulations during HIRI were revealed, highlighting the roles of Kupffer cells and monocyte-derived macrophages in modulating the hepatic microenvironment. The MIF and VISFATIN signaling pathways play important roles in the interaction between macrophages and other cells. These findings enhance our understanding of the mechanisms of PANoptosis in HIRI and provide a new basis and potential targets for prevention and treatment strategies for HIRI. Supplementary Information The online version contains supplementary material available at 10.1038/s41598-025-99264-6. Keywords: PANoptosis, Hepatic ischemia-reperfusion injury, Mononuclear phagocyte, Bioinformatics analysis, Machine learning Subject terms: Data mining, Functional clustering, Gene ontology, Gene regulatory networks, High-throughput screening, Machine learning, Microarrays, Predictive medicine, Virtual drug screening Introduction Hepatic ischemia-reperfusion injury (HIRI) is a pathophysiological process in which hepatocytes are damaged after ischemia, and the damage is aggravated after the restoration of blood flow^[38]1. It is commonly seen in liver transplantation, temporary blockage of blood flow during hepatic resection, trauma, and hemorrhagic shock^[39]2–[40]4. HIRI is characterized by a complex inflammatory response involving several key factors, including anaerobic metabolism, oxidative stress, mitochondrial dysfunction, intracellular calcium overload, Kupffer cells, neutrophils, cytokines, and chemokines^[41]5. Severe HIRI may lead to hepatocellular injury and necrosis or even liver failure, which may have adverse effects on the prognosis of patients. The mechanisms of HIRI have not been fully elucidated. The current methods used to prevent and treat HIRI, such as the use of antioxidants, the inhibition of inflammatory responses, and improvements in microcirculation, still lack satisfactory effects^[42]6. Therefore, understanding the pathogenesis mechanisms of HIRI is essential for developing effective preventive and therapeutic strategies. Previous studies have shown that the inflammatory response caused by HIRI can trigger the development of pyroptosis, apoptosis, necroptosis, ferroptosis, and other types of programmed cell death, which can exacerbate HIRI^[43]7–[44]9. PANoptosis is a newly identified form of programmed cell death that involves the simultaneous occurrence and regulation of pyroptosis, apoptosis, and necroptosis^[45]10,[46]11. Current research on the mechanisms related to PANoptosis has focused on neurologic IRI, whereas fewer studies have investigated PANoptosis in HIRI^[47]12–[48]14. Therefore, investigating the function of PANoptosis and its regulatory mechanisms in HIRI can deepen our understanding of the occurrence and development of HIRI and inspire the clinical application of HIRI prevention and treatment^[49]15. In this study, we comprehensively analyzed transcriptomic data before and after IRI during human liver transplantation. By combining bioinformatics analysis and machine learning methods, we identified PANoptosis-related biomarkers associated with HIRI risk and constructed and validated a predictive model for HIRI outcomes. In addition, the differences in the liver tissue microenvironment before and after IRI, especially the heterogeneity of the monocyte phagocyte system and intercellular communication, were explored. Results Analysis of PANoptosis-related differentially expressed genes Figure [50]1A illustrates that the ssGSEA scores for apoptosis, necroptosis, pyroptosis, and PANoptosis were significantly elevated in the HIRI group compared to controls. These findings suggest that these specific forms of programmed cell death might be activated during HIRI. The overlaps among the genes associated with apoptosis, necroptosis, and pyroptosis were shown in the Venn diagram (Fig. [51]S1A). Most of the GSEA results confirmed the activation of these pathways in the HIRI group (Fig. [52]S1B). PCA before and after batch correction demonstrated that our data integration effectively minimized batch effects (Fig. [53]1B). Performing DEG analysis between non-HIRI and HIRI groups revealed 2022 upregulated and 1313 downregulated genes, including 113 upregulated and 32 downregulated genes specifically associated with PANoptosis (Fig. [54]1C & Fig. [55]S1C). Gene enrichment analysis showed that DEGs were predominantly enriched in biological process-related GO terms such as regulation of protein kinase activity, mononuclear cell differentiation, leukocyte proliferation, and taxis. Additionally, KEGG pathway enrichment analysis indicated significant enrichment of DEGs in signaling pathways including PI3K-Akt, MAPK, TNF, IL-17, among others (Fig. [56]1D). Fig. 1. [57]Fig. 1 [58]Open in a new tab (A) Comparison of ssGSEA scores for PANoptosis-related genes between the non-HIRI and HIRI samples. (B) PCA of the GEO datasets before and after batch correction. (C) Volcano plot of DEGs between the non-HIRI and HIRI samples. Red indicates upregulated genes, and blue indicates downregulated genes. The first 5 upregulated and downregulated genes are labeled. (D) GO term and KEGG pathway enrichment analysis of 3335 DEGs between the non-HIRI and HIRI samples. Weighted gene coexpression network analysis of HIRI-related gene modules WGCNA was utilized to construct gene coexpression networks. Two outlier samples ([59]GSM3081915 and [60]GSM3081916) were identified and removed to ensure the quality of the network construction (Fig. [61]2A). Gene modules were identified based on the top 3,000 genes with a high median absolute deviation. When the soft threshold power was set to 6, the scale-free topology fit index reached 0.85, indicating that the network structure was optimal (Fig. [62]2B). Ultimately, six distinct gene modules represented by different colors were identified (Fig. [63]2C). Among these gene modules, the yellow gene modules presented the strongest correlation with the ssGSEA scores for the four programmed death signatures and HIRI outcomes (Fig. [64]2D). Additionally, a significant positive correlation was observed between the gene significance for HIRI and the module membership in the yellow module (Fig. [65]2E). Functional enrichment analysis of the yellow module revealed the involvement of the genes in a variety of biological processes and signaling pathways, including cytokine signaling in the immune system, vasculature development, VEGFA-VEGFR2 signaling, positive regulation of programmed cell death, and regulation of leukocyte activation (Fig. [66]2E). Fig. 2. [67]Fig. 2 [68]Open in a new tab (A) Two outliers are detected by sample clustering analysis. (B) Scale-free topological network analysis with various soft threshold powers in WGCNA. (C) Gene clustering dendrogram and corresponding 6 gene modules with different colors. (D) Correlation analysis between gene modules, ssGSEA scores of the programmed cell death gene signature, and HIRI outcome. (E) Scatterplot of correlation analysis between the importance of HIRI genes in the yellow module and module members. (F) Gene enrichment analysis of WGCNA identified HIRI-related genes. Identification and functional analysis of the hub genes To further refine the list of candidate genes, we took the intersection of 675 PANoptosis-related genes, 3335 DEGs, and 428 HIRI-related genes identified by WGCNA. The crossover identified a total of 44 common genes (Fig. [69]3A). These genes were subsequently analyzed via a protein-protein interaction network, and 23 genes with a minimum required interaction score of ≥ 0.9 were screened for subsequent study (Fig. [70]3B). The functional enrichment analysis revealed that these screened genes were mainly enriched in biological processes and pathways such as regulation of apoptotic signaling pathway, IL-17 signaling pathway, positive regulation of cytokine production, TNF signaling pathway, and response to unfolded proteins (Fig. [71]3C). Fig. 3. [72]Fig. 3 [73]Open in a new tab (A) Intersection of DEGs, PANoptosis-related genes, and HIRI-related candidate genes identified via WGCNA. (B) Protein-protein network analysis of HIRI-related gene intersections. (C) Gene enrichment analysis of HIRI-related gene intersections. (D) HIRI-related biomarkers identified by the LASSO algorithm. (E) HIRI-related biomarkers identified by the random forest algorithm. (F) HIRI-related biomarkers identified by the SVM-RFE algorithm. (G) HIRI-related biomarkers identified by the Boruta algorithm. (H) Intersection of HIRI-related biomarkers identified by four machine learning methods. (I) Differential expression of six HIRI-related genes between the non-HIRI and HIRI samples from the GEO combined dataset and [74]GSE151648. * P < 0.05, ** P < 0.01, *** P < 0.001. (J) ROC curves of the discrimination ability of six HIRI-related genes for HIRI outcome from the GEO combined dataset and [75]GSE151648. Key gene identification was performed using four machine-learning algorithms, including LASSO, RF, SVM-RFE, and Boruta. LASSO regression was used to select 13 candidate genes according to the minimum lambda value (Fig. [76]3D). The RF algorithm screened the top 10 genes according to the MeanDecreaseGini value ranking (Fig. [77]3E). The SVM-RFE algorithm identified 22 important genes (Fig. [78]3F), whereas the Boruta algorithm identified 16 important genes (Fig. [79]3G). The intersection of the genes identified by these algorithms was visualized via a Venn diagram, which revealed six common hub genes, namely, CEBPB, HSPA1A, HSPA1B, IRF1, SERPINE1, and TNFAIP3 (Fig. [80]3H). All six genes were highly expressed in the GEO combined dataset and [81]GSE151648 cohort (Fig. [82]3I), and they exhibited AUC values greater than 0.7 in predicting HIRI outcomes, indicating their potential as HIRI-related biomarkers (Fig. [83]3J). Establishment and performance evaluation of the nomogram A nomogram incorporating six hub genes for the prediction of HIRI was constructed and is presented in Fig. [84]4A. The ROC curve of this nomogram had an AUC value of 0.965, indicating high predictive accuracy. After performing 1000 bootstrap resamplings, the AUC values of the ROC curves remained consistently greater than 0.5, which was statistically significant (Fig. [85]4B). The calibration curve further confirmed the good agreement between the predicted probabilities and the actual observations (Fig. [86]4C). DCA demonstrated that the diagnostic model could provide greater net clinical benefit than the other diagnostic methods, such as the “test all” or “test none” approaches, over most risk thresholds ranging from 0.1 to 0.8 (Fig. [87]4D). This finding suggested the clinical value of the nomogram. In addition, the clinical impact curve illustrated that the nomogram could accurately predict HIRI outcomes as the risk threshold increased, further supporting its utility in clinical decision-making (Fig. [88]4E). Fig. 4. [89]Fig. 4 [90]Open in a new tab (A) Nomogram for HIRI prediction based on the transcriptome expression of six HIRI-related genes. (B) ROC curves of the nomogram after 1000 bootstrap resampling analyses. (C) Calibration curve of the nomogram after 1000 bootstrap resampling analyses. (D) Decision curve analysis of the nomogram. (E) Clinical impact curve analysis of the nomogram. (F) ROC curves of the nomogram in the [91]GSE151648 cohort. (G) The GEO combined dataset was divided into two clusters by consensus clustering analysis based on the six HIRI-related genes. (H) PCA of the clustering results for the GEO combined dataset. (I) Percentages of HIRI and non-HIRI samples in cluster 1 and cluster 2. The chi-square test results suggest that P < 0.01. (J) Comparison of ssGSEA scores for PANoptosis-related genes between the non-HIRI and HIRI samples in cluster 1 and cluster 2. (K) GSVA of cluster 1 and cluster 2. External validation of the diagnostic model in independent cohorts further confirmed its robustness and generalizability. The AUC value of the ROC curve was 0.9849 in the [92]GSE151648 cohort, which not only confirmed the high predictive performance of the model but also demonstrated that it was able to maintain good predictive ability across different datasets (Fig. [93]4F). Molecular typing based on the 6 HIRI-related hub genes The GEO combined dataset was classified into two clusters by consensus clustering analysis based on the expression patterns of the six HIRI-related hub genes. Cluster 1 contained 57 samples, while cluster 2 contained 95 samples (Fig. [94]4G). The results of the PCA demonstrated a significant clustering effect (Fig. [95]4H). Cluster 1 consisted of 55 non-HIRI samples and 2 HIRI samples, whereas cluster 2 consisted of 21 non-HIRI samples and 74 HIRI samples. Statistical analysis showed that the proportion of HIRI samples was significantly higher in cluster 2 than in cluster 1 (77.9% vs. 3.5%, P < 0.05) (Fig. [96]4I). The results of ssGSEA revealed that cluster 2 had significantly higher gene enrichment scores than cluster 1 for programmed cell death related to pyroptosis, apoptosis, necroptosis, and PANoptosis (P < 0.05) (Fig. [97]4J). GSVA further revealed the differences in biological processes between the two clusters. Cluster 1 was mainly enriched in NF-kappa B signaling, hypoxia, apoptosis, UV response, P53 pathway, inflammatory response, and many other biological processes. In contrast, cluster 2 was mainly enriched in fatty acid metabolism, protein secretion, bile acid metabolism, and oxidative phosphorylation (Fig. [98]4K). Function and interaction analysis of the 6 HIRI-related hub genes To further explore the functional roles of these six HIRI-related hub genes, the GSEA of the genes correlated with the expression of the hub genes was performed on the basis of the hallmark gene sets. The results revealed that the genes were involved in multiple complex biological processes, including TNFα signaling via NFκB, inflammatory response, allograft rejection, hypoxia, and apoptosis (Fig. [99]5A). Fig. 5. [100]Fig. 5 [101]Open in a new tab (A) Hallmark GSEA of the genes associated with the six HIRI-related genes. (B) UpSet plot of the potential molecular compounds for six HIRI-related genes. (C) Molecular docking simulation result of resveratrol and IRF1 (Affinity: − 6.1 kcal/mol). STRING constructed a protein interaction network for six HIRI-related genes, with IRF1 at its core. This positioning suggests that IRF1 may serve as a key regulatory molecule playing an essential role in the network (Fig. [102]S2A). Utilizing GOSemSim, we combined functionally annotated information from these six HIRI-related core genes within the GO database, encompassing molecular function, biological process, and cellular component categories. The geometric mean of the semantic similarity scores among these genes under three ontologies was calculated. The analysis results depicted in Fig. [103]S2B revealed that the proteins HSPA1A and transcription factor IRF1 exhibited the highest mean semantic similarity score, indicating their potential critical roles as core components within the six-gene interaction network. Furthermore, to predict the TF-miRNA-gene regulatory network of the hub genes, we utilized the RegNetwork database available on the NetworkAnalyst platform. In this regulatory network, CEBPB, IRF1, and TNFAIP3 emerged as central transcription factors, whereas HSPA1A, HSPA1B, and SERPINE1 were identified as core protein-coding genes. Notably, has-miR-155 displayed the most interactions with other genes within the network (Fig. [104]S2C). Small molecule compounds associated with the six genes in the context of reperfusion injury were searched through the CTD database. Cyclosporine, dexamethasone, estradiol, pirinixic acid, resveratrol, sodium arsenite, and tert-butylhydroperoxide were found to interact with six genes (Table [105]1; Fig. [106]5B). Resveratrol has been proven to attenuate HIRI in previous studies. Our molecular docking simulation results showed a potential interaction between resveratrol and IRF1 protein, providing new clues for further studies (Fig. [107]5C). Table 1. Potential molecular compounds associated with the six HIRI-related hub genes. Gene Inference network Inference score Reference count HSPA1A 4-phenylbutyric acid, Acetylcysteine, Antimycin A, Blood Glucose, cobaltous chloride, Cocaine, Cyclosporine, Deferoxamine, Dexamethasone, Edaravone, Estradiol, Ethanol, Hydrogen Peroxide, Indomethacin, Morphine, Oxygen, pirinixic acid, Quercetin, Ranitidine, Reactive Oxygen Species, Resveratrol, Rotenone, sodium arsenite, tert-Butylhydroperoxide, Troglitazone 55.59 79 HSPA1B Acetylcysteine, alpha-Tocopherol, Antimycin A, Blood Glucose, cinnamaldehyde, cobaltous chloride, Cocaine, Cyclosporine, Dexamethasone, Estradiol, Ethanol, Hydrogen Peroxide, Oxygen, pirinixic acid, Quercetin, Resveratrol, Rotenone, Sirolimus, sodium arsenite, tert-Butylhydroperoxide 46.51 66 IRF1 15-deoxy-delta(12,14)-prostaglandin J2, Acetylcysteine, Antimycin A, Cyclosporine, Dexamethasone, Estradiol, Ethanol, Flutamide, Glucose, Hydrogen Sulfide, mangiferin, pirinixic acid, Reactive Oxygen Species, Resveratrol, Rosiglitazone, Simvastatin, sodium arsenite, Streptozocin, tert-Butylhydroperoxide, Troglitazone 51.71 51 TNFAIP3 Acetylcysteine, Adenosine, Antimycin A, bardoxolone methyl, Cyclosporine, Dexamethasone, Dihydrotestosterone, Estradiol, Ethanol, Gallic Acid, Hydrogen Peroxide, Mannitol, Oxygen, pirinixic acid, Quercetin, Resveratrol, Rotenone, sodium arsenite, tert-Butylhydroperoxide, Testosterone, Troglitazone 49.37 63 CEBPB 4-phenylbutyric acid, Acetylcysteine, Antimycin A, baicalein, Capsaicin, cobaltous chloride, Cyclosporine, Dexamethasone, Epinephrine, Estradiol, Ethanol, Gallic Acid, Glucose, Hydrogen Peroxide, Indomethacin, Mifepristone, Nicotine, Nitroprusside, Oxygen, pirinixic acid, Quercetin, Reactive Oxygen Species, Resveratrol, Rosiglitazone, Rotenone, Sevoflurane, sodium arsenite, Streptozocin, Telmisartan, tert-Butylhydroperoxide, Testosterone, Troglitazone, Vitamin D 78.54 84 SERPINE1 15-deoxy-delta(12,14)-prostaglandin J2, Blood Glucose, Capsaicin, cinnamaldehyde, cobaltous chloride, Cyclosporine, Dexamethasone, Dihydrotestosterone, Epinephrine, Estradiol, Fenofibrate, Flutamide, Glucose, Heparin, Hydrogen Peroxide, Hydrogen Sulfide, Indomethacin, Losartan, mangiferin, Mifepristone, Morphine, Nicotine, notoginsenoside R1, N-(oxo-5,6-dihydrophenanthridin-2-yl)-N, N-dimethylacetamide hydrochloride, Oxygen, Pentoxifylline, pirinixic acid, Quercetin, Ranitidine, Resveratrol, Rosiglitazone, Rotenone, Simvastatin, Sirolimus, sodium arsenite, Streptozocin, Superoxides, Tacrolimus, Telmisartan, tert-Butylhydroperoxide, Testosterone, Troglitazone, Vitamin D 98.93 93 [108]Open in a new tab Annotation and functional analysis of cell subpopulations The scRNA-seq data from [109]GSE171539 were processed and analyzed via the standard Seruat workflow. Cluster analysis with a resolution of 0.8 resulted in the identification of 18 cell clusters, which were subsequently annotated into nine distinct cell subpopulations, including mononuclear phagocytes, NK cells, T cells, B cells, endothelial cells, plasma cells, hepatocytes, fibroblasts, and cholangiocytes (Fig. [110]6A). Figure [111]6B illustrates the expression of canonical marker genes across these different cell subclusters, providing clear evidence for their identity. The heatmap in Fig. [112]6C displays the expression levels of the top five DEGs in each cell subpopulation, highlighting the unique transcriptional profiles of these cell types. Fig. 6. [113]Fig. 6 [114]Open in a new tab (A) Cell subpopulation clustering analysis and cell type annotation of the [115]GSE171539 scRNA-seq data. (B) Classic marker gene expression in different cell populations. (C) Heatmap of the expression of the top 5 genes in different cell populations. (D) Heatmap of different cell populations in the non-HIRI and HIRI samples. (E) Heatmap of the proportions of different cell populations in the non-HIRI and HIRI samples. (F) Metabolism pathway activity in the non-HIRI and HIRI samples. (G) Metabolism pathway activity in different cell populations. (H) Comparison of AUCell scores for PANoptosis-associated gene sets in the non-HIRI and HIRI samples. * P < 0.05, ** P < 0.01, *** P < 0.001. (I) Comparison of AUCell scores for PANoptosis-associated gene sets in the non-HIRI and HIRI samples from different cell populations. * P < 0.05, ** P < 0.01, *** P < 0.001. By comparing the proportions of each cell subpopulation between the HIRI and non-HIRI groups, we found that the proportion of mononuclear phagocytes was significantly increased in the HIRI group, which may reflect the enhanced inflammatory response caused by mononuclear phagocytes and other cells (Fig. [116]6D and E). Metabolic pathway analysis revealed that metabolic pathways such as oxidative phosphorylation, inositol phosphate metabolism, glycolysis/gluconeogenesis, and glutathione metabolism were more active in the HIRI group, whereas citrate cycle (TCA cycle) was more active in the non-HIRI group (Fig. [117]6F). In addition, hepatocytes were identified as the most metabolically active cell subpopulation during HIRI (Fig. [118]6G). The results of the GO term biological processes and KEGG pathways enrichment analysis of the nine cell subpopulations are shown in Fig. [119]S3, indicating that these cell subpopulations are involved in a wide range of biological processes, including the immune response, energy metabolism, cell survival, and tissue repair after HIRI. The AUCell scores for the programmed cell death gene sets were calculated and compared between the HIRI and non-HIRI groups to further explore the activity of PANoptosis in different cell subpopulations. The results showed that the AUC scores of the four types of programmed cell death were significantly higher in the HIRI group than in the non-HIRI group, indicating increased activity of these biological processes (Fig. [120]6H). Specifically, in the HIRI group, all the cell subpopulations presented increased AUCell scores indicating PANoptosis, suggesting that these cells were more actively involved in PANoptosis. (Fig. [121]6I). Mononuclear phagocyte reclassification and exploration of the cell trajectory Figure [122]7A shows that the expression levels of the 6 HIRI-related core genes were significantly greater in the liver tissues of the HIRI group than in those of the non-HIRI group. Notably, similar results were observed in mononuclear phagocytes (Fig. [123]7B). Mononuclear phagocytes were reclustered for further analysis. Unsupervised clustering at a resolution of 0.2 identified two distinct macrophage subpopulations (Fig. [124]7C). These two subclusters were annotated based on the gene expression patterns of their highly upregulated signature genes. Cluster 0 (monocyte-derived macrophages, MoMFs) highly expressed S100A8, S100A9, and S100A12, similar to the gene expression profiles of circulating monocytes in peripheral blood. Cluster 1 (Kupffer cells, KCs), on the other hand, expressed tissue-resident cell marker genes, such as C1QA, C1QB, and C1QC, suggesting that they were Kupffer cells in the liver (Fig. [125]7D). Fig. 7. [126]Fig. 7 [127]Open in a new tab (A) Distribution and expression of six HIRI-related genes in the non-HIRI and HIRI samples. * P < 0.05, ** P < 0.01, *** P < 0.001. (B) Differential expression of six HIRI-related genes in mononuclear phagocytes from the non-HIRI and HIRI samples. * <0.05, ** <0.01, *** <0.001. (C) Reclustering analysis and cell type annotation of mononuclear phagocytes in the non-HIRI and HIRI samples. (D) Heatmap of the top 10 genes expressed in KCs and MoMFs. (E) The proportion of MoMFs and KCs in the non-HIRI and HIRI samples. The chi-square test results were statistically significant (P < 0.05). (F) Comparison of AUCell scores for PANoptosis-associated gene sets in the non-HIRI and HIRI samples * P < 0.05, ** P < 0.01, *** P < 0.001. (G) Pseudotemporal trajectories of mononuclear phagocytes. (H) Expression levels of six HIRI-related genes along the pseudotemporal trajectory of mononuclear phagocytes. (I) Heatmap of the top 100 DEGs along the pseudotemporal trajectory. A comparison of the proportions of different macrophage subpopulations revealed a significant increase in the proportion of MoMFs in the HIRI group compared with that in the non-HIRI group, suggesting greater influx of MOMFs during HIRI (Fig. [128]7E). GO term enrichment analysis of the marker genes for MOMF revealed that these genes were primarily involved in biological processes such as regulation of actin filament organization, response to type II interferon, plasminogen activation, and regulation of actin cytoskeleton organization. However, the marker genes of the KC subpopulation were enriched in biological processes such as leukocyte-mediated immunity, antigen processing and presentation of exogenous peptide antigens, regulation of T-cell activation, regulation of leukocyte cell-cell adhesion, and lymphocyte-mediated immunity. According to the KEGG pathway enrichment analysis, the marker genes of MOMFs were mainly involved in leukocyte transendothelial migration, amino acid biosynthesis, carbon metabolism, glycolysis and gluconeogenesis, and shigellosis, whereas the marker genes of KCs were associated with phagosomes, lysosomes, antigen processing, and presentation, type I diabetes mellitus, and allograft rejection (Fig. [129]S4A). A comparison of the PANoptosis AUCell scores between the two macrophage subpopulations revealed that both subpopulations presented increased PANoptosis activity in the HIRI groups (Fig. [130]7F). Figure [131]7G shows the pseudotime change trajectory of the macrophage subpopulation transition from MOMFs to KCs. As inferred from the CytoTRACE2 results, the pseudotemporal trajectory should be from MoMF to KC (Figure [132]S4B). It could be hypothesized that the proportion of MoMFs converted to KCs might gradually increase as HIRI progresses. As shown in Fig. [133]7H, the gene expression levels of HSPA1A and HSPA1B gradually increased with pseudotime, whereas the expression of other genes, especially IRF1 and TNFAIP3, decreased and then increased. These findings provide valuable insights into the cellular and molecular mechanisms underlying HIRI and highlight potential targets for therapeutic intervention. Figure [134]7I demonstrates the changes in DEGs along the pseudotemporal trajectory, where HSPA1A and HSPA1B were highly expressed in cluster 2. GO term enrichment analysis showed that these DEGs were mainly enriched in the biological processes, including leukocyte chemotaxis and migration. KEGG pathway enrichment analysis revealed that these DEGs were predominantly enriched in IL-17, TNF, and NF-κB signaling pathways, as well as antigen processing and presentation (Fig. [135]S4C). Intercellular communication between mononuclear phagocytes and other cells CellChat was used to analyze intercellular communication in the scRNA-seq data. The total number of cell subpopulations was lower in the HIRI group than in the non-HIRI group. However, the overall strength of intercellular communication was significantly increased in the HIRI group (Fig. [136]S5A). Figure [137]8A presents heatmaps of the inferred differences in the number and strength of interactions among all the cell subpopulations. Specifically, KCs, MoMFs, NK cells, T cells, and endothelial cells exhibited significantly greater intercellular communication in the HIRI group. Fig. 8. [138]Fig. 8 [139]Open in a new tab (A) Differences in the number and strength of intercellular communication between the non-HIRI and HIRI groups. The red color represents an increase in the HIRI group compared with the non-HIRI group, and the blue color represents a decrease in the HIRI group compared with the non-HIRI group. (B) Comparison of different signaling pathways involved in intercellular communication in the non-HIRI and HIRI groups. Red represents the signaling pathways enriched in the non-HIRI group, and blue represents the signaling pathways enriched in the HIRI group. Stacked bar plots are shown on the left, and nonstacked bar plots are shown on the right. (C) Signaling pathways accompanied by the upregulation of ligand-receptor gene expression in the HIRI group. (D) Visualization of the VISFATIN signaling pathway in the non-HIRI and HIRI groups. Comparison of signal senders and receivers of the VISFATIN signaling pathway between the non-HIRI and HIRI groups. (E) Visualization of the MIF signaling pathway in the non-HIRI and HIRI groups. Comparison of signal senders and receivers of the MIF signaling pathway between the non-HIRI and HIRI groups. The scatter plots highlighted substantial changes in the intensity of afferent and efferent signaling pathways for certain cell subpopulations, with MoMFs and KCs showing particularly notable changes (Fig. [140]S5B). Fig. [141]S5C compares and visualizes the signaling pathways with large variations in afferent and efferent signals for the intercellular communication of KCs and MOMFs, respectively. The VISFATIN and MIF signaling pathways strongly differ in outgoing interaction strength between mononuclear phagocytes and other cells. Figure [142]8B shows comparative bar plots of the conserved and specific signaling pathways in the non-HIRI and HIRI groups. Several signaling pathways, including VISFATIN, LIFR, IL4, EGF, IL6, PDGF, IL1, ANNEXIN, and FASLG, were enriched in the HIRI group. Additionally, the GH, SEMA3, RESISTIN, SAA, and ANGPTL signaling pathways were uniquely expressed in the HIRI group. Figure [143]8C presents the signaling pathways associated with upregulated gene expression of ligand-receptor pairs between KCs, MOMFs, and other cells in the HIRI group compared with those in the non-HIRI group. These results suggested that a complex and dynamic intercellular communication network during HIRI, especially the VISFATIN and MIF signaling pathways (MIF-ACKR3, MIF-(CD74 + CXCR4), NAMPT-(ITGA5 + ITGB1), and NAMPT-INSR), was upregulated in the HIRI group. The VISFATIN and MIF signaling pathways were notably enhanced in macrophages and other cell subpopulations in the HIRI group compared with those in the non-HIRI group. In addition, the role of macrophages as important signal senders of the VISFATIN and MIF signaling pathways was more prominent in the HIRI group than in the non-HIRI group (Fig. [144]8D and E). These findings further highlight the critical role of macrophages in the HIRI microenvironment and provide inspiration for exploring HIRI-associated macrophage therapeutic targets. Discussions HIRI is a major cause of graft dysfunction following liver transplantation and can affect the recovery of patients and even lead to death^[145]16. At the initiation of HIRI, ischemia and hypoxia lead to disruption of the mitochondrial electron transport chain and the generation of reactive oxygen species (ROS), mitochondrial permeability transition, and intracellular calcium overload, leading to hepatocyte and endothelial cell injury. Damaged cells can release damage-associated molecular patterns to activate pattern recognition receptors in the immune system. Subsequently, pro-inflammatory cytokines and chemokines were produced by KCs to recruit inflammatory cells, such as monocytes and neutrophils, to the injury site and exacerbate the inflammatory response. Finally, ROS and other toxic metabolites can be generated during reperfusion to further damage the cell membrane structure and function, leading to cell death or the induction of programmed cell death, such as apoptosis and ferroptosis^[146]17,[147]18. PANoptosis, a newly discovered form of programmed cell death, including pyroptosis, apoptosis, and necroptosis, has been proven to be associated with IRI in multiple organs, including heart, lungs, and brain^[148]19–[149]21. Although pyroptosis, apoptosis, and necroptosis are involved in HIRI, the occurrence and regulation of PANoptosis in HIRI have been less well studied^[150]22. Intriguingly, the consistent results of bulk RNA-seq and scRNA-seq bioinformatics analyses suggested that PANoptosis-related genes are enriched in HIRI. In this study, six genes associated with PANoptosis were identified by comprehensively analyzing the transcriptomic data of HIRI patients before and after liver transplantation. The expression patterns and functional roles of these genes could provide new insights into their regulatory mechanisms in HIRI. Gene set enrichment analysis revealed that the genes related to these genes were regulated mainly by NF-κB. The NF-κB pathway is closely related to the pathology of HIRI. In addition, NF-κB also interacts with PI3K/Akt, MAPK, and Nrf2 signaling pathways to regulate HIRI. Inhibition of the NF-κB pathway can delay the inflammatory response and cell death. Therefore, the NF-κB pathway may be a potential target for HIRI^[151]23. Resveratrol has been shown in vitro cellular experiments to modulate genes involved in NF-κB signaling, including IRF1^[152]24. Corresponding molecular docking simulation was provided in our study. IRF1 and TNFAIP3 are the transcription factors associated with HIRI^[153]25. IRF1 is a member of the interferon regulatory factor family. IRF1 expression is mainly induced by IFN and promoted by the activation of the iNOS/NO signaling pathway during HIRI. Studies have shown that IRF1 not only promotes hepatocyte apoptosis and necrosis during HIRI but also further exacerbates liver injury by inducing the activation and chemotaxis of various immune cells^[154]26. TNFAIP3 is a negative regulator of the NF-κB pathway, which helps attenuate liver injury induced by inhibiting inflammatory responses^[155]27. However, other studies have suggested that the overexpression of TNFAIP3 may exacerbate HIRI^[156]28. The transcriptional activator CCAAT/enhancer binding protein β, encoded by CEBPB, is important for regulating immune and inflammatory responses and may significantly influence the progression of HIRI. Studies have shown that CEBPB expression can be induced by proinflammatory cytokines such as IL-1, IL-6, IL-8, and TNF. CEBPB can promote the production of proinflammatory cytokines and acute-phase proteins to exacerbate liver injury^[157]29,[158]30. In addition, an association between C/EBPβ and IRI in experimental animal models has been demonstrated. HSF1 is activated and upregulates the expression of heat shock proteins (HSPs) under stress conditions. HSPA1A and HSPA1B encode Hsp70-1 and Hsp70-2, respectively. They play a key role in maintaining cellular homeostasis and protecting cells from injury. In a renal IRI model, the upregulation of HSPs was closely associated with injury and recovery of renal function^[159]31. The molecular chaperone function and antiapoptotic, antioxidant, and anti-inflammatory effects of Hsp70 are important for attenuating liver injury. SERPINE is a key serine protease inhibitor that regulates fibrinolysis. It acts primarily by inhibiting tissue-type plasminogen activator and urokinase-type plasminogen activator^[160]32. Meanwhile, HSF1 was able to upregulate SERPINE expression in endothelial cells under stress conditions^[161]33. SERPINE can promote the infiltration of neutrophils into damaged areas and enhance the inflammatory response^[162]34. This can lead to further tissue damage. In our study, miR-155 was identified as the key miRNA regulator of six HIRI-related genes. It has been shown that miR-155 has a regulatory role in liver diseases, including HIRI, and that the inhibition of miR-155 can alleviate HIRI^[163]35. Our findings increase the understanding of PANoptosis-related molecular mechanisms in HIRI and provide potential therapeutic targets. A novel diagnostic model was constructed based on of the six PANoptosis-related core genes, and its high predictive performance in the internal and external validation cohorts was confirmed. In addition, we delved into the characteristics of molecular typing based on these six core genes to provide an important basis for individualized treatment. In recent years, the development of single-cell sequencing technologies has greatly facilitated the understanding of the complex tissue microenvironment. In the context of HIRI, changes in the immune microenvironment are particularly critical. In this study, we revealed the heterogeneity and dynamics of macrophage subpopulations during HIRI via scRNA-seq analysis. Macrophages constitute the largest population of intrinsic immune cells in the liver and are highly heterogeneous, playing an important role in maintaining liver homeostasis and regulating disease processes. In our study, mononuclear phagocytes were divided into KCs and MOMFs according to the previous single-cell transcriptome map^[164]36. KCs originate from the progenitor cells of the embryonic yolk sac, reside in liver tissues, and maintain self-renewal. When the liver is damaged, KCs can secrete inflammatory cytokines and chemokines to promote immune cell infiltration into the liver. Among them, monocytes derived from the bone marrow are recruited to the liver and derived into MoMFs^[165]37. Pseudotemporal trajectory analysis indicated that MOMFs are gradually transformed into KCs. Studies have shown that MoMFs can differentiate into KCs after the depletion of KCs^[166]38. This process is accompanied by increased expression of heat shock proteins, suggesting the resolution of liver injury. Further analysis revealed that intercellular communication via MIF-ACKR3, MIF-(CD74 + CXCR4), NAMPT-(ITGA5 + ITGB1), and NAMPT-INSR ligand-receptor pairs in macrophages may have a significant effect on HIRI. The upregulation of NAMPT-(ITGA5 + ITGB1) and NAMPT-INSR between macrophages and endothelial cells was also observed in a chronic hepatitis B (CHB) mouse model. These findings suggest that macrophage modulation of endothelial cell dysfunction and angiogenesis is present in both HIRI and CHB^[167]39. MIF is a proinflammatory cytokine that promotes inflammation, apoptosis, oxidative stress, and angiogenesis. MIF binds to the receptors CD74/CD44, CXCR2, CXCR4, and CXCR7, and activates downstream MAPK, AMPK, and Akt signaling to promote inflammatory responses^[168]40. Our findings not only reveal the dynamic changes in macrophage subpopulations and their functional properties during HIRI but also provide important clues for further research on the pathological mechanisms of HIRI and the development of new therapeutic strategies. However, this study has several limitations. First, the sample size of this study was relatively small because there are fewer publicly available datasets on HIRI. However, a detailed search of the GEO database was conducted to ensure that sufficient data were available to meet sample size requirements and that all the samples had complete information. Second, our study was conducted via bioinformatics analysis, and cellular and animal experiments were not performed. However, we combined bulk RNA-seq and scRNA-seq data analysis with multiple machine-learning methods to identify biomarkers and construct prediction models. Finally, future experiments should be conducted to validate these important findings. In conclusion, in this study, we identified six genes associated with PANoptosis and constructed a high-performance HIRI prediction model by comprehensively analyzing transcriptomic data before and after IRI in human liver transplantation. In addition, we utilized single-cell sequencing technology to explore the heterogeneity and dynamic changes in macrophage subpopulations during HIRI in detail. These findings not only deepen our understanding of the mechanism of PAN apoptosis in HIRI but also provide a new theoretical basis and potential targets for future clinical applications. Future studies can further explore the value of these biomarkers and signaling pathways in clinical practice to provide more effective strategies for the prevention and treatment of HIRI. Materials and methods Data collection To identify essential biomarkers associated with IRI during liver transplantation, this study comprehensively analyzed bulk RNA sequencing (bulk RNA-seq) data from human liver tissues before and after IRI from the Gene Expression Omnibus (GEO) database. A dataset consisting of seven datasets ([169]GSE87487, [170]GSE113024, [171]GSE112713, [172]GSE12720, [173]GSE14951, [174]GSE15480, and [175]GSE23649) was constructed as the development cohort, and [176]GSE151648 was used as an independent external validation cohort to confirm the findings. The R package sva was used to address batch effects and conduct data integration. Principal component analysis (PCA) was used to assess batch removal. A total of 675 PANoptosis-related genes were obtained after combining the genes associated with apoptosis, pyroptosis, and necroptosis^[177]41,[178]42. These genes were selected to explore the underlying molecular mechanisms of IRI. In addition, to explore the cellular heterogeneity in the HIRI microenvironment in greater depth, the single-cell RNA sequencing (scRNA-seq) dataset [179]GSE171539 was included in this study. The detailed data sources and clinical information for the above datasets are summarized in Table [180]2. Table 2. Data sources for bulk RNA-sequencing, single-cell RNA-sequencing transcriptomic data, and matched clinical information for the non-HIRI and HIRI samples. GEO accession Platform Non-HIRI sample (pre-liver transplantation, n) HIRI sample (post-liver transplantation, n) Reperfusion time (h) [181]GSE151648 [182]GPL21290 23 23 2 [183]GSE87487 [184]GPL11154 4 4 2 [185]GSE113024 [186]GPL10999 3 3 1 [187]GSE112713 [188]GPL14951 11 11 1 [189]GSE12720 [190]GPL570 21 21 1 [191]GSE14951 [192]GPL570 5 5 2–3 [193]GSE15480 [194]GPL6244 6 6 1.5 [195]GSE23649 [196]GPL6947 33 33 2 [197]GSE171539 [198]GPL20795 1 1 2 [199]Open in a new tab Single-sample gene set enrichment analysis To explore the differences in the four types of programmed cell death (PANoptosis, apoptosis, pyroptosis, and necroptosis) between non-HIRI and HIRI states, we calculated and compared the enrichment scores of each sample in the GEO combined dataset via the single-sample gene set enrichment analysis (ssGSEA) algorithm. Differentially expressed gene analysis Differentially expressed gene (DEG) analysis between 76 non-HIRI samples and 76 HIRI samples from the GEO combined dataset was performed via the limma R package. The thresholds of significant difference were set as an absolute value of log fold change (logFC) > 0 and a false discovery rate (FDR) < 0.05. Weighted gene coexpression network analysis Weighted gene coexpression network analysis (WGCNA) was used to construct gene coexpression networks and identify gene modules highly correlated with PANoptosis ssGSEA scores and HIRI outcomes. This gene module identified by WGCNA intersected with the PANoptosis-related DEGs. The shared genes were uploaded to the STRING platform ([200]https://cn.string-db.org/) for protein interaction analysis and to the Metascape platform ([201]https://metascape.org/) for gene functional annotation analysis. Identification of HIRI-related biomarkers Four machine learning algorithms, including least absolute shrinkage and selection operator (LASSO) regression, random forest (RF), support vector machine recursive feature elimination (SVM-RFE), and Boruta, were used for diagnostic biomarker selection in the GEO merged dataset. The expression levels of the selected HIRI-related biomarkers were compared between non-HIRI and HIRI samples. The predictive performance of the biomarkers for HIRI outcome was evaluated by calculating metrics such as the area under the receiver operating characteristic curve (ROC). Establishment and performance evaluation of the prediction model A nomogram based on the transcriptome expression levels of HIRI-related biomarkers was constructed via logistic regression and validated in internal and external validation cohorts for its ability to predict HIRI outcomes. The predictive performance of the nomogram was evaluated by area under the ROC curve (AUC), calibration curve, decision curve analysis (DCA), and clinical impact curve analyses. Consensus clustering analysis Consensus clustering analysis of the GEO combined dataset was conducted on the basis of the expression patterns of the screened HIRI-related genes. The clustering effect was evaluated via PCA. The percentages of non-HIRI and HIRI samples in different clusters were compared via the chi-square test. The ssGSEA scores for the four types of programmed cell death in different clusters were compared. Gene set variation analysis (GSVA) was applied to estimate and compare differences in hallmark pathway activity between different clusters. Functional analysis of HIRI-related biomarkers Activated hallmark pathways associated with HIRI-related gene expression were identified via gene set enrichment analysis (GSEA). Gene interaction networks and transcription factor (TF)-microRNA (miRNA) coregulatory networks of the selected HIRI-related biomarkers were visualized via STRING ([202]https://cn.string-db.org/) and NetworkAnalyst ([203]https://www.networkanalyst.ca/). The semantic similarities of the GO terms associated with the biomarkers were calculated and compared using GOSemSim R package. Prediction of potential molecular compounds The Comparative Toxicogenomics Database (CTD) ([204]https://ctdbase.org/) was used to confirm the HIRI-related biomarkers further and predict their corresponding potential molecular compounds. The structure data files of the molecular compounds were downloaded from the PubChem database ([205]https://pubchem.ncbi.nlm.nih.gov/). The protein information of the core biomarkers was obtained from GeneCards database ([206]https://www.genecards.org/) and UniProt database ([207]https://www.uniprot.org/). The protein sequences and 3D structures were downloaded from AlphaFold protein structure database ([208]https://alphafold.ebi.ac.uk/). Molecular docking was simulated with AutoDockTools 1.5.7 and visualized with PyMOL 3.0. Single-cell transcriptomic data analysis Cellular heterogeneity and unique functional states were revealed through the Seurat workflow. Genes expressed in fewer than three cells were excluded from quality control. Meanwhile, cells expressing fewer than 500 genes with unique molecular identifier counts of less than 1000 were excluded. In addition, cells whose content of mitochondrial genes was greater than 25%, the content of ribosomal genes greater than 50%, and the content of erythroid genes greater than 1% were excluded from the study. Normalization was performed using SCTransform. Dimensionality reduction was performed using PCA. Data integration and batch effect correction were performed by the Harmony algorithm. Distinct cell clusters were identified by uniform manifold approximation and projection (UMAP) clustering analysis and annotated based on the results of the Single R package and differential gene expression profiles. The AUCell scores of the four types of programmed death were estimated and compared between the HIRI and HIRI groups. The differential expression of the screened HIRI-related core genes was verified in the scRNA-seq data. In addition, metabolic pathway activity comparison and Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analyses were performed on different cellular subpopulations. Monocle 2 R package was used to infer the dynamic cell states and cell transitions in mononuclear phagocytes during HIRI^[209]43. CytoTRACE2 was used to predict the cellular differentiation potential of monocyte macrophage subpopulations, validating the results of pseudotemporal trajectory analysis^[210]44. CellChat R package was employed to explore and compare potential differences in intercellular communication between the non-HIRI and HIRI groups, and relevant signaling pathways were visualized, especially in mononuclear phagocytes^[211]45. Electronic supplementary material Below is the link to the electronic supplementary material. [212]Supplementary Material 1.^ (34.4MB, tif) [213]Supplementary Material 2.^ (2MB, tif) [214]Supplementary Material 3.^ (35.9MB, tif) [215]Supplementary Material 4.^ (35.6MB, tif) [216]Supplementary Material 5.^ (33.2MB, tif) [217]Supplementary Material 6.^ (16.1KB, docx) Author contributions Z.C. contributed to the data analysis and drafted the manuscript. J.Z. and L.Z. contributed to data verification. Y.L. and T.Z. contributed to data acquisition and preprocessing. Z.C., J.Z., L.Z., Y.L., T.Z., and X.S. contributed to the data analysis. Y.X. and X.L. contributed to the study design and revision of the manuscript. All the authors read and approved the final manuscript. Funding This work was supported by National High Level Hospital Clinical Research Funding (2022-PUMCH-C-049, 2022-PUMCH-A-237). Data availability Data is provided within the manuscript or supplementary information files. Declarations Competing interests The authors declare no competing interests. Footnotes Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Contributor Information Yiyao Xu, Email: xuyiyao@pumch.cn. Xin Lu, Email: luxin@pumch.cn. References