Identification of Renal Ischemia-Reperfusion Injury Subtypes and Predictive Model for Graft Loss after Kidney Transplantation Based on Programmed Cell Death-Related Genes

Clinically, ischemia-reperfusion injury (IRI) is an unavoidable complication following renal transplantation, leading to acute kidney injury, delayed graft function (DGF), and failure of the transplanted kidney [1‒3]. Substantial clinical research has demonstrated that IRI is merging as a pivotal contributor to deterioration in long-term graft survival, negatively affecting the survival rate and the quality of patients’ lives [4, 5]. However, due to the complexity and diversity of IRI pathophysiology, as well as the high heterogeneity of patient clinical symptoms and long-term prognosis, prediction of patients’ prognosis remains difficult and challengeable. Thus, molecular subtypes and systematical identification of molecules related to prognosis may help solve the huge challenge of patient heterogeneity of renal IRI and further facilitate to uncover novel targeted therapies for kidney injury.

IRI is an intricate pathophysiological process, initiating a myriad of cellular responses aimed at either repairing or leading to injury [6, 7]. The failure of these repair mechanisms, resulting in cell death, emerges as the most critical endpoint and the aspect most intimately linked to prognostic outcomes, hence garnering widespread attention. With deeper insights into the molecular mechanisms of cell death, programmed cell death (PCD) was known as a regulated cell death mediated by dedicated molecular machines, which plays pivotal roles in both physiological homeostasis and disease pathology [8, 9]. To date, eighteen distinct modalities of PCD have been identified, encompassing anoikis, alkaliptosis, apoptosis, autophagy, cuproptosis, entotic cell death, entosis, immunogenic cell death, ferroptosis, lysosome-dependent cell death, methuosis, necroptosis, netotic cell death, NETosis, oxeiptosis, paraptosis, parthanatos, and pyroptosis [10‒12]. Different patterns of cell death indicate different biological outcomes. IRI can cause PCD by releasing endogenous molecules termed damage-associate molecular patterns (DAMPs) to trigger a series of inflammatory and cytotoxic responses, which finally results in an aseptic inflammatory reaction [13]. Several PCD patterns such as necroptosis, apoptosis, and autophagy have been reported that are positively associated with poor prognosis after renal transplantation [14]. However, there are no relevant reports on the relationship between 18 PCD patterns and IRI after kidney transplantation. Therefore, a comprehensive understanding of the PCD patterns after IRI is conductive to deepen the understanding of transplanted kidneys, explore the determinants or influencing factors of prognosis and develop an optimal predictive model for prognosis.

In the current study, we performed the first research to comprehensively analyze the relationship between 18 PCD patterns and IRI, propose cell death types that affect prognosis and establish a predictive model for graft loss after kidney transplantation. In order to solve these problems, we first classified IRI patients into four categories by consensus clustering and explored different signatures of cell death, metabolism, immune infiltration and clinical characteristics in these subtypes. Next, we identified the statistically significant patterns of cell death that were predictive to prognosis after renal transplantation by various machine learning methods and screened out the differentially expressed genes (DEGs) in PCD which were associated with graft loss. Finally, we developed a robust model for the prediction of graft loss after renal transplantation based on 39 hub genes from 18 PCD patterns, which may help with risk stratification, early detection, and timely intervention of patients with poor prognosis after kidney transplantation.

Figure 1 illustrates the step-by-step procedures of our study. First of all, we downloaded eight datasets of biopsied kidney samples related to IRI after kidney transplantation from the Gene Expression Omnibus (GEO) database [15]. Eighteen patterns of PCD and the key regulatory genes were extracted from previous reports [16]. Consensus clustering was then conducted to stratify IRI patients on the basis of PCD-related genes to determine distinct PCD-related IRI subtypes. Subsequently, differential enrichment analysis was applied to analyze the differences in cell death patterns, metabolic signatures, and immune characteristics among the different PCD subtypes, and the prognostic differences indicated by graft function loss were assessed using the Kaplan-Meier (K-M) survival analysis. Next, statistically significant PCD patterns that correlated to poor prognosis were screened by the intersection of results from three machine learning algorithms. The DEGs significantly related to graft function loss in each PCD pattern were defined and a predictive model of graft loss was constructed by PCD index (PCDI). The K-M curve and time-dependent receiver operating characteristic (ROC) curves were used to validate the sensitivity and specificity of our predictive model.

The gene expression data and the clinical information of patients were obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) [15]. Altogether eight datasets of biopsied transplanted kidney samples were chosen for the subsequent analysis and validation, including GSE21374 (n = 282), GSE43974 (n = 554), GSE37838 (n = 78), GSE10419 (n = 30), GSE53769 (n = 36), GSE90861 (n = 46), GSE61739 (n = 96), and GSE30718 (n = 47). A detailed description of these datasets is provided in online supplementary Table 1 (for all online suppl. material, see https://doi.org/10.1159/000540158). Data filtering, background correction, quality control, and log2 transformation of these datasets were performed. Normalization was then conducted by “normalizeBetweenArrays” function in the “limma” R package. Additionally, we merged the datasets and applied a batch correction using the Combat method from the “sva” R package.

Based on literature research, we gathered key regulatory genes of 18 PCD patterns [17]. A total of 1,964 PCD-related genes including 7 alkaliptosis genes, 338 anoikis genes, 580 apoptosis genes, 367 autophagy genes, 19 cuproptosis genes, 15 entotic cell death genes, 23 entosis genes, 87 ferroptosis genes, 34 immunogenic cell death genes, 220 lysosome-dependent cell death genes, 8 methuosis genes, 101 necroptosis genes, 8 netotic cell death genes, 24 NETosis genes, 5 oxeiptosis genes, 66 paraptosis genes, 9 parthanatos genes, and 52 pyroptosis genes were gathered. After removing 416 duplicates, 1,548 PCD-related genes were finally used for our analysis.

We conducted a consensus clustering analysis to identify IRI subtypes based on PCD-related genes (IRI PCD subtypes). Merged expression matrix was a subset to gather the 1,548 PCD-related cluster genes for our analysis. The transplanted kidney samples post-IRI were then classified into distinct subtypes by employing the “Consensus ClusterPlus” package [18]. To determine the optimal number of consensus clustering, we performed a test by running the algorithm by setting k from 2 to 10 sequentially. The determination of the optimal k was dependent on the best consensus among the items in each cluster. In addition, we employed principal component analysis (PCA) to visualize the clustering patterns of the samples.

Gene Set Variation Analysis (GSVA) represents a gene set enrichment analysis through an unsupervised and nonparametric approach that estimates the score attributed to a particular pathway or signature based on transcriptomic data [19]. We acquired 18 PCD patterns and 84 metabolism signatures from previously published studies [16, 17]. By setting the method to single sample gene set enrichment analysis (ssGSEA) in the “GSVA” R package, we calculated gene set scores for each sample corresponding to these 18 PCD patterns and 84 metabolism signatures.

To assess the degree of immunological infiltration comprehensively, various bioinformatic algorithms including ESTIMATE [20], CIBERSORT [21], Xcell [22], MCPcounter [22], Quantiseq [23], and EPIC [24] were employed. Each algorithm employed a unique strategy and gene expression signature to estimate the abundance of different immune cell subpopulations. The ESTIMATE R package was used to quantify the relative abundance of immune and stromal cells between the IRI PCD subtypes based on their gene expression profiles. The CIBERSORT, Xcell, MCPCounter, Quantiseq, and EPIC algorithms were employed to calculate the relative abundance indicated by the infiltration proportion of each immune cell type. Finally, we utilized the Kruskal-Wallis test to determine the significance of differences in immune infiltration. Heatmap and boxplot analyses were used to provide a visual representation of immune cell proportions in the kidney allografts with significant enrichment.

Three machine learning methods including PCA [25], ssGSEA [26], and z-score [27] were utilized to calculate the gene set score of each PCD pattern by “IOBR” R package [28] with a selection criterion of p value <0.05, and the PCD patterns significantly correlated with prognosis were screened and further defined by Univariate Cox regression analysis. The “VennDiagram” package was used to explore the overlapping cell death patterns, respectively, defined in the three machine learning algorithms.

First, prognosis-associated genes in 18 PCD patterns were screened using univariate Cox regression (hazard ratio (HR) ≠ 1 and p value <0.05). Subsequently, the candidate genes were further generated from least absolute shrinkage and selection operator (LASSO) regression feature selection and multivariate Cox regression model using the R package “glmnet.” For visualization, forest plots of selected genes were plotted using R package “forestplot.” Finally, we established the risk model and divided all the kidney samples into high- and low-risk groups based on the risk score. The risk scores of each PCD model were calculated using a linear combination of the regression coefficient and the expression of candidate genes (Riskscore = $∑i=1ngeneCoefi×geneExpi$⁠). A detailed formula for each risk model was addressed in online supplementary Table 3. Different expressions of the screened genes were compared between the high- and low-risk groups.

In order to construct a prognostic model with high efficacy, 101 unique combinations from ten diverse machine learning algorithms [29] were integrated to develop PCDI. These machine learning algorithms included LASSO, CoxBoost, Elastic Net, Random Survival Forest (RSF), Stepwise Cox, Support Vector Machine (SVM), Ridge, Super Partial Correlation (SuperPC), Gradient Boosting Machine (GBM), and Partial Least Squares with Cox regression (plsRcox). PCDI was generated by the combination of RSF and Ridge. The model’s robustness was validated by the single PCD pattern set and four external sets (Bi et al. [30], Dou et al. [31], Wu et al. [32], and Chen et al. [33]). Patients from GSE21374 dataset were categorized into high-risk and low-risk groups according to the score of PCDI, and the prognostic differences between the two groups were defined using K-M curves. Additionally, ROC curves were generated to evaluate the PCDI’s prognostic efficacy.

Statistical analyses and scientific graphing were performed by R (version 4.3.2). The distributions of continuous and nominal variables were characterized by descriptive statistics. Student’s t test and Mann-Whitney U nonparametric test were applied to normal and abnormal distribution variables, respectively. Between-group comparisons were conducted via the Wilcoxon test. To evaluate the impact of different cell death patterns on survival outcomes, Cox regression models and K-M survival analysis were utilized. ROC analysis was performed using the “pROC” package, and the area under the curve was calculated. A significance level of p < 0.05 was considered statistically significant.

In order to define the phenotypes of renal PCD following ischemia-reperfusion kidney injury, we took advantage of the available eight GEO datasets which include 755 biopsied transplanted kidney tissues and conducted cluster analysis based on the expression of 1,548 PCD-related genes. The corresponding clinical information of the eight GEO datasets is summarized in online supplementary Table 1. As shown in Figure 2a and b, after removing the batch effect, the consensus clustering analysis classified the gene expression profiles of the transplanted kidney samples into distinct subtypes. The cumulative distribution function plot displayed the fluctuation trend of k = 2–10 (Fig. 2c). Based on the cophenetic, dispersion and silhouette metrics, k = 4 was determined as the optimal number of clusters, which had the minimum fluctuation in cumulative distribution function curve and clear boundaries in consensus matrix heatmap (Fig. 2c, d). Therefore, the transplanted kidney samples post-IRI were ultimately divided into four PCD subtypes (IRI PCD subtypes) and were designated as PCD-A (n = 245), PCD-B (n = 167), PCD-C (n = 125), and PCD-D (n = 218), respectively (Fig. 2e). Detailed sample counts information for each subtype was provided in online supplementary Table 2. The PCA indicated a remarkable difference in the whole transcriptomics gene expression among the four IRI PCD subtypes (Fig. 2e). Gene expression differences were observed by heatmap analysis (Fig. 2f).

We next investigated the clinical relevance of the four IRI PCD subtypes in the eight GEO datasets. As shown in Figure 2g, the four IRI PCD subtypes could be detected in nearly all the datasets, except for an absence of PCD-B in GSE30718 (with a relatively limited sample size, n = 47). In GSE43974 and GSE10419 datasets where the donor types had been informed, we found that in the kidneys from donors after brain death or donors after cardiac death, PCD-A and PCD-B subtypes were the main cell death patterns, while in the kidneys from living donors, PCD-D was the predominant cell death pattern (Fig. 2h). Additionally, we attempted to explore the relationship between the four IRI PCD subtypes and kidney prognosis. In the GSE21374 dataset which includes a prospective cohort with long-term follow-up, we found that patients classified as PCD-D had the best kidney prognosis among the four subtypes (p < 0.01), while those defined as PCD-A exhibited the worst long-term graft survival rate compared with patients with the other three IRI PCD subtypes (Fig. 2i).

We next investigated the features of cell death in the four IRI PCD subtypes. The gene set scores of the 18 PCD patterns were quantified by GSVA and visualized in the clustered heatmap (Fig. 3a). Boxplot was painted to quantify the gene set score of subtype-specific cell death patterns by differential enrichment analysis (Fig. 3b). As shown in Figure 3, each IRI PCD subtype was composed of various cell death patterns. Specifically, PCD-A exhibited enrichment in more cell death signatures compared with other subtypes, including the upregulated expression of NETosis, immunogenic cell death, anoikis, necroptosis, entotic cell death, and ferroptosis (Fig. 3a). PCD-B was observed to have higher scores of pyroptosis, apoptosis, autophagy, alkiliptosis and methuosis. Signatures of netotic cell death, lysosome-dependent cell death, paraptosis, entosis, and oxeiptosis showed stronger enrichment in PCD-C (Fig. 3a). Notably, PCD-D exhibited an extensive downregulation of cell death except entotic cell death, cuproptosis, and parthanatos, which was in line with its relevance to a better prognosis (Fig. 2i). Kruskal-Wallis test was used to compare statistical differences, and the statistical results displayed significant differences of these four subtypes among the 18 PCD patterns (Fig. 3b).

Given that the kidney is an organ with a highly active metabolism, we further investigated whether the different IRI PCD subtypes exhibited varying metabolic signatures. Altogether 84 metabolism processes from nine metabolic modules were measured by GSVA. As shown in the heatmap of Figure 4a, the metabolic signature analysis revealed a gradually upregulated expression of metabolism pathways from PCD-A to PCD-D. GSVA score was used to visualize subtype-specific metabolic signatures as previously mentioned (Fig. 4b). Notably, the metabolism signatures of glycan biosynthesis and metabolism were the most distinguishing features of PCD-C and metabolism of other amino acids was highly enriched in PCD-D (Fig. 4a). PCD-C and PCD-D showed a higher score of seven metabolic modules, including carbohydrate metabolism, energy metabolism, lipid metabolism, nucleotide metabolism, amino acid metabolism, metabolism of cofactors and vitamins, and metabolism of terpenoids and polyketides, which revealed an extensive activation of metabolic pathways in these two subtypes (Fig. 4a). Additionally, PCD-A exhibited a general differentially downregulated pattern of multiple metabolic modules, accompanied by a worse outcome after kidney transplantation (Fig. 2i, 4a).

To elucidate the overall differences in immune characteristics among the four IRI PCD subtypes, the ESTIMATE algorithm was applied to calculate the immune and stromal scores. The immune and stromal score showed a remarkably higher level in PCD-A than the other three subtypes (p < 0.001), whereas PCD-D exhibited the lowest immune and stromal score across the four subtypes (p < 0.001; Fig. 5a, b). We further observed the expression of immune checkpoints. Compared with other subtypes, PCD-C displayed a significantly higher expression of several immune checkpoint genes, including TLR9, CTLA4, and PDCD1. The expression of CD274 was significantly higher in PCD-A (Fig. 5a, c). These immune checkpoints with significant variation may serve as targets for immunotherapy. Additionally, immune infiltration was examined to depict the immunological landscape based on five different immune-related databases through GSVA (Fig. 5a). We quantified the abundance of multiple immune microenvironment cells in boxplots (Fig. 5d–h). In detail, PCD-A presented significantly highest abundance of various immune cell populations (including activated NK cells, monocytes, M1 macrophage, mast cell, aDC, cDC, cytotoxic lymphocytes, myeloid dendritic cells, B cells, M2 macrophage, NK cells, CD8⁺ T cells, and CD4⁺ T cells, etc.) compared with other three subtypes (p < 0.001). The lowest immune infiltration was investigated in PCD-D, mainly containing naïve B cells, activated CD4⁺ memory T cells, and monocyte lineage (Fig. 5d–h). Of note, PCD-A demonstrated a significantly higher abundance of fibroblasts, which seemed to be consistent with the poor prognosis in this subtype analyzed above (Fig. 5f, h).

Three different machine learning algorithms were conducted to screen the specific PCD patterns related to prognosis, which was indicated by the loss of graft function. Using rigorous statistical analysis and cut-off values, we further narrowed down the cell death patterns closely related to prognosis. As depicted in Figure 6a–c, 14 cell death patterns were defined as significantly correlated with graft loss by the PCA algorithm (Fig. 6a), 13 patterns were identified through ssGSEA algorithm (Fig. 6b), and 16 patterns were recognized by the z-score calculation (Fig. 6c). The overlapping 11 cell death patterns (immunogenic cell death, necroptosis, pyroptosis, entosis, entotic cell death, apoptosis, anoikis, paraptosis, cuproptosis, methuosis, and NETosis) of the above three algorithms with statistical significance (p < 0.05) were shown by the Venn diagram (Fig. 6d). In the GSE21374 dataset with long-term follow-up, K-M curves were plotted for these 11 types of cell death to evaluate their relevance to prognosis (Fig. 6e, f). The risk score was calculated for each case of kidney sample, and all the cases were divided into high-risk and low-risk groups according to the median risk score. Based on Hazard Ratio values, immunogenic cell death was the optimal pattern indicating a poor prognosis. Additionally, cell death patterns such as necroptosis, pyroptosis, entosis, entotic cell death, apoptosis, and anoikis also showed a positive association with graft function loss (Fig. 6f).

We further screened the key regulatory genes of the 18 cell death patterns that were significantly associated with graft function loss. As displayed by the forest plots in Figure 7a, all of the immunogenic cell death-related genes (BAX, CD8A, EIF2AK3, HMGB1, IFNG, IL1B, IL6, and NT5E) and the entotic cell death-related genes (CTNNA1, CYBB) were positive predictors of graft function loss, while both oxeiptosis-related genes including KEAP1 and AIFM1 and parthanatos-related genes including MIF and OGG1 were negatively correlated with graft function loss. The other 14 cell death patterns contained both positively and negatively correlated genes with graft loss. Specifically, anoikis-related genes including ITGA6, ATF4, PIK3R3, MUC1, EEF1A1, RHOB, apoptosis-related genes including ITM2C, POLB, TGFB1, TNFRSF12A, lysosome-dependent cell death-related genes including CTSO, DNASE2, GGA3, HEXB, STXBP1, necroptosis-related genes including CAPN2 and STAT6, and pyroptosis-related genes including CHMP4C, IL1A, NOD1, PLCG1, TNF were mainly positively correlated with graft loss. On the contrary, autophagy-related genes including CPTP, KLHL3, ORMDL3, ferroptosis-related genes including NQO1, PCBP2, PEBP1, PGD, TFRC, and methuosis-related genes including GIT1 and MTOR were negatively correlated with graft loss (Fig. 7a). Based on these key regulatory genes, multivariate Cox regression was utilized to establish risk models of each cell death pattern, and GSE21374 dataset with long-term follow-up was used to validate the efficiency of the risk models. Patients were divided into high- and low-risk groups according to the best cut-off values of risk score. All differential expression of key regulatory genes in the 18 cell death patterns between high-risk and low-risk groups was shown in Figure 7b.

We next aimed to establish a risk model for predicting long-term graft loss based on the cell death-related hub genes. By intersecting DEGs of graft function loss and significant genes inferred from Cox proportional hazards regression, a set of 39 genes derived from the 18 cell death patterns were identified as risk genes linked to graft function loss and were considered to be suitable for prognostic model development (Fig. 8a). Forest plots were utilized to visualize the significant prognostic linkage of these genes to graft loss as indicated by the hazard ratio values (Fig. 8b). In order to construct a prognostic model with high efficacy, 101 unique combinations from ten distinct machine learning algorithms were utilized. Among them, the combination of RSF and Ridge algorithms exhibited optimal performance, and the risk score of this model was then named PCDI (Fig. 8c). In the GSE21374 dataset with long-term follow-up, PCDI demonstrated the best efficacy when compared with the four published prediction models (Bi et al. [30], Dou et al. [31], Wu et al. [32], and Chen et al. [33]) and any risk model constructed by each PCD pattern (Fig. 8d). Patients were divided into high-risk and low-risk groups according to the score of PCDI. The K-M curve manifested that recipients in the high-risk group had significantly worse long-term graft survival (p < 0.0001, Fig. 8e). Time-dependent ROC curves showed the area under the curve values of 0.5-, 1-, 2-, 3-, and 4-year graft survival were 0.888, 0.91, 0.926, 0.923, and 0.923, respectively, indicating that our predictive model had good sensitivity and specificity in predicting graft loss after kidney transplantation (Fig. 8f). The risk score of each recipient was calculated to divide patients into high- and low-risk type (Fig. 8g). Recipients with higher risk score were more probably to undergo graft failure (Fig. 8h). Among the 39 genes mentioned above, 4 genes (NFS1, ATP6V0A4, ATP6V0D2, and ATP6V1C2) were upregulated in the low-risk group and 35 genes were elevated in the high-risk group (Fig. 8i).

IRI poses a significant risk to kidney transplants, potentially leading to graft loss and transplant failure [34, 35]. There is mounting evidence that the IRI process activates various PCD pathways in transplanted kidneys [14, 36]. Nevertheless, due to the huge heterogeneity in the IRI pathophysiology during the complicated clinical process of kidney transplantation, the cell death patterns in the kidney tissues remain unclear. It is imperative to explore cell death patterns that determine prognosis in order to enhance therapeutic strategies and improve long-term outcomes. This study is the first to thoroughly examine the association between PCD and IRI, along with its impact on the prognosis of transplanted kidneys. Through consensus clustering, four distinct IRI PCD subtypes with different features of cell death, metabolism, and immune infiltration were identified. Distinct relevance of IRI PCD subtypes to the long-term prognosis was defined, and the key regulatory genes of variant cell death were explored as predictors of graft function loss. A predictive model was constructed with the highest sensitivity and specificity among existing models, offering a promising tool for forecasting graft survival following kidney transplantation.

As our understanding deepens, an increasing body of evidence highlights the role of cell death in IRI after kidney transplantation, especially for necroptosis [37] and ferroptosis [38]. The emergence of these cell death patterns often indicated negatively impacting both short- and long-term post-transplant outcomes [39]. Following IRI, DAMPs were released from necrotic cells and subsequently activated inflammasomes, thereby driving tissue inflammation at the initial stage of necroptosis [40]. Lau et al. [41] investigated that receptor-interacting protein kinase 3 (RIPK3)-mediated necroptosis exacerbated donor kidney inflammatory injury and reduced allograft survival in mouse models. Chen et al. [33] explored the occurrence of ferroptosis during IRI after kidney transplantation both in cellular and mouse models, which eventually led to severe inflammatory infiltration and worse outcomes. However, the available studies mainly focused on the impact of individual PCD forms on renal IRI. Systematic analysis of the relationship between all reported PCD patterns and IRI, thereby exploring the specific cell death pattern that determines prognosis remains blank. Cluster analysis serves as a powerful tool to define subtypes of IRI in kidney transplant patients, based on similar cell death signatures. Therefore, our research has a significant impact on identifying patients’ subtypes and assisting in precise therapy.

Among the four IRI PCD subtypes defined in our study, PCD-A is associated with the poorest kidney prognosis. This subtype exhibits significant enrichment in various PCD patterns and is characterized by a pronounced cell death signature related to the immune response, notably including NETosis, immunogenic cell death, and necroptosis. It is reported that various cell death patterns can be strong and potentially lethal triggers of innate immune system activation and vice versa [40]. NETosis involves the expulsion of neutrophil extracellular traps into the extracellular matrix to combat pathogens directly, playing a vital role in the innate immune response [42]. Immunogenic cell death is a type of PCD that elicits an immune reaction through the release of DAMPs from dying cells, thereby attracting immune cells to the site of cell death [43]. Necroptosis-mediated inflammation occurs via several mechanisms including release of cytokines and DAMPs, RIP3- and mixed lineage kinase domain-like protein (MLKL)-independent activation of inflammasome [44]. Notably, interaction between cell death and inflammation can give rise to an auto-inflammation loop, leading to exaggerated cell death and further inflammation – a process referred to as necroinflammation. This could trigger systemic inflammation and potentially result in graft injury and rejection post-transplantation [40]. Consistently, our findings also defined immunogenic cell death, with immune response as the main feature, as the optimal form to predict poor prognosis. Other immune response-related cell death patterns, such as necroptosis, pyroptosis, apoptosis, and NETosis, have likewise been strongly associated with adverse outcomes. Nevertheless, in contrast to the other three subtypes, PCD-D is characterized by extensive downregulation of cell death and the lowest level of immune infiltration, along with a more favorable prognosis. This suggests that PCD-D might represent a distinct subtype of PCD. Studies have confirmed that kidney transplants from living donors, which typically involve shorter ischemia times, lead to better long-term graft survival and are therefore considered a preferable option for transplantation [45]. Correspondingly, the primary donor source for PCD-D is living donors, which correlates with a milder ischemic process and less activation of cell death pathways. Therefore, blocking cell death-inflammation loop could be crucial in developing specific therapeutic interventions.

Recent clinical studies have indicated that metabolic failure was the primary effector mechanism of renal IRI [46]. Our study suggested the distinct metabolic characteristics across different IRI PCD subtypes. Notably, PCD-A characterized by extensive cell death patterns and the poorest prognosis exhibited significant metabolic paralysis with a broad downregulation of various metabolic modules, including carbohydrate and energy metabolism. On the contrary, a general activation of catabolic pathways was observed in PCD-D, corresponding to the best prognosis among the four subtypes. Metabolic dysfunction is known to correlate with graft function deterioration after kidney transplantation [47]. Lindeman et al. [46] described immediate and persistent post-reperfusion metabolic paralysis, including normoxic glycolysis and persistent post-reperfusion ATP/GTP catabolism, as key drivers of IRI that potentially lead to DGF. Tricarboxylic acid (TCA) cycle defect at the level of the oxoglutarate-dehydrogenase complex with selective post-reperfusion release of α-ketoglutarate and impaired recovery of tissue succinate levels was identified in clinical DGF cases [48]. Moreover, metal-dependent cell deaths, such as ferroptosis, have been implicated in exacerbating cellular damage through their involvement in various metabolic processes including carbon transformation, protein synthesis, the TCA cycle, and glycolysis [49‒51]. Collectively, there is an interactive relationship and reciprocal causation between metabolism condition and cell death. The in-depth study of the metabolic characteristics of different IRI PCD subtypes is of great significance for understanding the mechanism by which IRI reduced long-term renal graft function.

Our study defined 42 genes positively correlated with prognosis and 27 genes negatively correlated with prognosis among 18 PCD patterns. Most of these genes have been reported as key regulatory genes involved in cell death. For example, DAMPs released during immunogenic cell death such as ATP, high-mobility group box 1 (HMGB1), and heat shock proteins (HSPs) can promote antigen presentation and immune activation [12]. Catenin alpha 1 (CTNNA1) has been found to facilitate the process of entosis [52]. It has been reported that NOD1 plays fundamental and pleiotropic roles in host defense against infection and the control of inflammation, becoming upregulated during autophagy and pyroptosis [53, 54]. Additionally, Tang et al. [13] screened three ferroptosis-related hub genes (IL6, ATF3, and JUN) that were closely relevant to immune response and might be diagnostic biomarkers and therapeutic targets for IRI. Coincidentally, these genes mentioned above were associated with poor prognosis in our study. Hence, identifying changes in the expression of these genes may indicate the activation or suppression of corresponding cell death patterns, offering new insights for targeted therapy in clinical practice.

To effectively predict long-term graft function loss after kidney transplantation, we next constructed a prognostic model based on 39 PCD hub genes. Compared with any single model established by each cell death pattern and all four published prediction models (Bi et al. [30], Dou et al. [31], Wu et al. [32], and Chen et al. [33]), our predictive model demonstrated the best efficacy to date. There is some crosstalk among distinct PCD phenotypes, which makes it difficult to obtain satisfactory results by blocking a single type of PCD [55, 56]. Thus, our predictive model provides a significant advantage for kidney transplant patients in future clinical applications. First, IRI patients with different cell death profiles can be classified into different subtypes using our model, facilitating stratified diagnosis and guiding precise treatment. Second, based on different cell death forms, patients can be divided into high-risk and low-risk groups, allowing for the identification of high-risk patients with poor prognosis and timely intervention to improve outcomes. Finally, for patients predicted to have poor prognosis by our model, we may identify specific PCD genes contributing to poor outcomes. This enables us to define the critical cell death patterns in the transplanted kidney, potentially serving as therapeutic targets for individualized precision therapy.

There are some limitations in our study. First, the kidney prognosis data are sourced exclusively from a single dataset (GSE21374), which might introduce a selection bias. Therefore, the efficiency and stability of our prediction model require prospective validation across broader cohorts. Moreover, while we have identified many key regulatory genes associated with graft loss post-kidney transplantation, further experimental studies are still needed to validate these findings. Lastly, considering the significant protein degradation associated with cell death, future research should ideally incorporate proteomics to further validate and analyze the cell death patterns across different subtypes, enabling more precise conclusions.

In conclusion, this study is the first to systematically describe and define distinct PCD subtypes in IRI patients after kidney transplantation by employing data from various sources and synthesizing clinically consistent cohorts. Our findings broaden the existing knowledge of the correlation between graft loss and different cell death patterns in these patients. The prediction model, which focuses on specific PCD-related genes, holds strong potential for predicting graft loss, potentially aiding in risk stratification for patients post-transplantation. Future studies could focus on targeting the specific cell death pattern to develop precise treatment strategies for transplanted kidney.

An ethics statement was not required for this study type since no human or animal subjects or materials were used.

The authors have no conflicts of interest to declare.

This study was supported by grants from the National Natural Science Foundation of China (Nos. 82130021, 8220031230), Michigan Medicine-PKUHSC Joint Institute for Translational and Clinical Research (BMU2022JI002), the Beijing Young Scientist Program (BJJWZYJH01201910001006), Capital’s Funds for Health Improvement and Research (CFH2022-1-4071, 2020-JKCS009), CAMS Innovation Fund for Medical Sciences (2019-I2M-5-046).

Jing Ji, Yuan Ma, and Xintong Liu collected, analyzed, and interpreted the data. Jing Ji and Zehua Li drafted the manuscript. Xintong Liu, Qingqing Zhou, Xizi Zheng, and Ying Chen revised the manuscript. Zehua Li and Li Yang conceived the study, interpreted the data, and revised the manuscript. Li Yang supervised the study. All authors approved the final manuscript.

Jing Ji, Yuan Ma and Xintong Liu contributed equally to this work.

All data generated or analyzed during this study are included in this article and its online supplementary material. Further inquiries can be directed to the corresponding author.