Integrated Analysis and Validation of Ferroptosis-Related Genes Associated with Ischemia/Reperfusion Injury in Lung Transplantation

Introduction

Lung transplantation is the only therapeutic option to save patients with end stage lung disease.¹ However, ischemia/reperfusion injury of the graft during this process is one of the main causes of primary graft dysfunction, which can further increase the risk of chronic allograft dysfunction and mortality.² Therefore, alleviating lung ischemia/reperfusion injury is crucial for improving the prognosis of lung transplant recipients.

The term of ferroptosis was coined in 2012,³ a novel cell death form driven by iron-dependent phospholipid peroxidation, distinguishing it from other cell death types such as apoptosis, necrosis, pyroptosis, and autophagy.⁴ The essence of ferroptosis lies in the depletion of glutathione, with a decrease in the activity of glutathione peroxidase 4 (GPX4). Lipid peroxides cannot be metabolized via the GPX4-catalyzed glutathione reductase reaction, leading to the oxidation of lipid by ferrous iron ions and subsequent generation of reactive oxygen species, thereby triggering ferroptosis.⁴ In recent years, increasing research has found that ferroptosis occurs during organ ischemia/reperfusion injury and is recognized as a significant mode of cell death in this process.^5–10 Zhao et al observed activation of ferroptosis-related signaling pathways during ischemia/reperfusion in lung transplantation using human lung biopsy tissues.¹¹ Furthermore, a finding demonstrated that liproxstatin-1, a ferroptosis inhibitor, could alleviate lung ischemia/reperfusion injury. Additionally, studies have identified beta-aminoisobutyric acid,¹² salidroside,¹³ hydroxycitric acid,¹⁴ and irisin¹⁵ can ameliorate lung ischemia/reperfusion injury by inhibiting ferroptosis. Therefore, inhibiting ferroptosis may serve as a target for ameliorating lung ischemia/reperfusion injury in lung transplantation. In ischemia/reperfusion injury of the lungs, comprising ischemia and reperfusion phases, ferroptosis has been shown to occur in both processes through in vitro experiments utilizing BEAS-2B cells.¹¹ However, the exact mechanism of ferroptosis in lung ischemia/reperfusion injury remains largely unclear. Hence, further investigation into the precise mechanism is crucial, as it can provide valuable insights for identifying potential biomarkers and therapeutic targets related to ferroptosis in lung ischemia/reperfusion injury.

High-throughput sequencing technologies provide a powerful approach for studying the pathogenic characteristics of ferroptosis in lung ischemia/reperfusion injury. In this study, we obtained pre- and post-reperfusion data from lung transplant patients through the GEO database. By intersecting with ferroptosis-related genes, we identified FRGs, as well as module genes identified by WGCNA algorithm. Two machine learning algorithms, least absolute shrinkage and selection operator (LASSO) and support vector machine-recursive feature elimination (SVM-RFE) were utilized to analyze overlapping genes, thereby uncovering potential key genes associated with ferroptosis in lung ischemia/reperfusion injury. Furthermore, we established a murine model of lung ischemia/reperfusion injury for validation. Mendelian randomization is an epidemiological analysis method that utilizes genetic variation to determine the causal relationship between exposure factors and outcomes. In this study, we further elucidated the causal role of candidate genes in the occurrence of PGD using a two-sample Mendelian randomization approach. Finally, consensus clustering analysis based on FRGs identified distinct subtypes, offering valuable insights for early prognostic evaluation and the development of personalized treatment strategies for lung transplant patients.

Materials and Methods

Data Resources

The two datasets used in this study, GSE145989 and GSE18995, were sourced from the GEO database (http://www.ncbi.nlm.nih.gov/geo, accessed on February 28th, 2024). GSE145989 comprises 67 pairs of matched pre-reperfusion and post-reperfusion human lung transplant samples, while GSE18995 consists of 18 pre-reperfusion and 17 post-reperfusion samples (Table S1). Since both datasets consist of microarray data, sample size calculations were conducted using the “pwr” package in R, a method commonly applied in microarray studies to ensure statistical robustness. Based on an effect size (Cohen’s d = 0.5), a significance level of 0.05, and a target statistical power of 0.8, the required sample size was estimated to be at least 63 samples per group. These results confirmed that the sample size was sufficient for robust analysis. Additionally, 564 ferroptosis-related genes were obtained from the FerrDb database (https://www.zhounan.org/ferrdb/current/, accessed on February 28th, 2024).

Analysis of Differentially Expressed Genes (DEGs)

We employed the R package “limma” to identify DEGs between pre-reperfusion and post-reperfusion transplanted lungs in the GSE145989 dataset. The statistical criteria for differential expression were an |log2-fold change (logFC)| > 0.585 and adjusted P-value < 0.05. Visualization of downregulated and upregulated DEGs was performed using the R package “ggplot2” for volcano plots. The intersection of DEGs with ferroptosis-related genes named FRGs was illustrated using the Venn diagram (http://www.ehbio.com/test/venn/#/). Visualization of the expression levels of FRGs in pre-reperfusion and post-reperfusion transplanted lungs was conducted using the R package “pheatmap” for heatmap generation.

Functional Enrichment Analysis

To analyze the biological functions and pathways associated with FRGs, the R package “clusterProfiler” was used to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses. In the GO analysis, identified FRGs were categorized into three classes: biological processes (BP), cellular components (CC), and molecular functions (MF). KEGG enrichment analysis was conducted to predict signaling pathways. The criterion for significant enrichment was P < 0.05. The “ggplot2” package in R software was utilized to generate scatter plots.

Wgcna

The R package “WGCNA” was employed to identify module genes associated with reperfusion of lung transplantation. A Venn diagram (http://www.ehbio.com/test/venn/#/) was utilized to illustrate the shared genes between FRGs and module genes.

Identification of Key Genes

Further screening was conducted to identify genes associated with ferroptosis in lung ischemia/reperfusion injury using LASSO logistic regression. The R package “glmnet” was employed for LASSO logistic regression, where the minimum lambda was considered optimal. Through the SVM-RFE algorithm with the R packages “e1071” and “caret”, we identified key genes by pinpointing the point of lowest cross-validation error. Subsequently, receiver operating characteristic (ROC) curves and area under the curve (AUC) calculations were performed on the GSE18995 dataset using the R package “pROC” to evaluate the sensitivity and specificity of key genes.

The Mouse Model of Lung Ischemia/Reperfusion Injury

8-to 12-week-old male C57BL/6 mice were purchased from the Animal Center of Renmin Hospital of Wuhan University. The experimental protocol was approved by the Ethics Committee of Renmin Hospital of Wuhan University (permit NO.20221203). In order to establish lung ischemia/reperfusion injury (LIRI) model, mice were anesthetized with isoflurane and mechanically ventilated using an animal ventilator (PhysioSuite^TMRoVent, KENT, USA). The ventilator settings were adjusted to a respiratory rate of 150 breaths/min, a tidal volume of 0.4 mL, and a peak inspiratory pressure of less than 20 cmH₂O. Heparin (20U/kg, Solarbio, China, Cat# H8060) was administered via external jugular vein injection. A left thoracotomy was performed by cutting the left fourth rib to expose the left hilum. The left pulmonary hilum was clamped using a microvascular clip for 1 hour, followed by release of the clamp to initiate reperfusion. After 4 hours, mice were euthanized with an overdose of pentobarbital sodium, and the left lung was harvested for subsequent analysis. During the surgical intervention, rectal temperature was maintained between 36.5°C and 37.5°C. Sham mice underwent only thoracotomy without pulmonary hilum occlusion. The ferroptosis inhibitor ferrostatin-1 (Fer-1, 5 mg/kg, MedChemExpress, USA, Cat# HY-100579) was administered via intraperitoneal injection at the onset of reperfusion.

Histological Analysis

The left lung was placed in 4% paraformaldehyde at room temperature for 24 hours to fix the tissue. After fixation, the tissue was dehydrated through a gradient ethanol series, cleared in xylene, and then embedded in paraffin. The paraffin-embedded lung tissue was sectioned at 5 μm thickness using a microtome. The sections were then mounted on glass slides for Hematoxylin and Eosin (HE) staining. In the HE staining procedure, the sections were first stained with hematoxylin, followed by eosin staining to highlight the cytoplasm and extracellular matrix. After staining, the sections were dehydrated again, mounted with neutral balsam, and photographed. Pathological scoring of lung injury was carried out according to reference 16.¹⁶

Real-Time Quantitative Polymerase Chain Reaction (RT-qPCR)

Total RNA was extracted from mouse lung tissues using Trizol reagent (Servicebio, China, Cat# G3013). After RNA extraction, 1 μg of total RNA was reverse transcribed into cDNA using a cDNA synthesis kit (Vazyme, China, Cat# R223-01). The reverse transcription reaction was performed according to the manufacturer’s instructions. RT-qPCR was conducted using Power Up SYBR Green Master Mix (Servicebio, China, Cat# G3326-05) on a Roche LightCycler (Roche, Switzerland). The reaction mixture consisted of 10 μL SYBR Green Master Mix, 1 μL of cDNA, 0.4 μM forward and reverse primers, and RNase-free water to a total volume of 20 μL. The amplification conditions included an initial denaturation at 95°C for 5 min, followed by 40 cycles of 95°C for 15s, 60°C for 30s, and 72°C for 30s. The relative expression levels of target genes were normalized to β-actin, and the fold changes were calculated using the 2^−ΔΔCt method. The primer sequences utilized are listed below (Table 1).

Table 1 Gene-Specific Primers Used for RT-qPCR

Iron Detection

Lung tissues were homogenized in PBS buffer (0.1 M, pH 7.4) at a 1:9 (weight) ratio and centrifuged at 3500 rpm for 10 minutes at 4°C to remove debris. The iron concentration was determined using an iron assay kit (Nanjing Jiancheng Bioengineering Institute, China, Cat# A039-2-1) and measured at 520 nm using a microplate reader (Thermo Fisher Scientific, USA).

Malondialdehyde (MDA) and Reduced Glutathione (GSH) Assays

The levels of MDA and GSH in lung tissues were determined using assay kits (Nanjing Jiancheng Bioengineering Institute, China, Cat# A003-1-2 and A006-1-1), following the manufacturer’s protocols. The MDA content was measured at 532nm and the GSH content at 405nm, both using a microplate reader (Thermo Fisher Scientific, USA).

A Predictive Model and Nomogram of LRI

A predictive model for LIRI was developed using logistic regression. The “rms” was employed to create a nomogram, while the “rmda” package was used to generate decision curve analysis (DCA) and calibration curves.

Immune Cell Infiltration Analysis

The CIBERSORT algorithm and single-sample gene set enrichment analysis (ssGSEA) was utilized to infer the infiltration status of immune cell types in the samples. “Heatmaps” and “box plots” in R were employed to visualize the differences in immune cell infiltration between pre-reperfusion and post-reperfusion samples, with P < 0.05 considered statistically significant. Additionally, Pearson correlation analysis was conducted to validate the association between key genes and immune infiltration. Visualization of the results was performed using the R package “corrplot”.

MR Analysis

The MR analysis aims to investigate the causal relationship between candidate genes and the risk of PGD, where the instrumental variables are single nucleotide polymorphisms (SNPs). The cis-eQTL data for genes were obtained from eQTLGen, and the genome-wide association study (GWAS) summary statistics for the outcome event was sourced from the Integrated Epidemiology Unit (IEU) database (Table S2). The cis-eQTL data, with a p-value threshold of 5×10⁻⁶, were used as exposure factors for the MR analysis. If no relevant SNPs were available, the corresponding cis-eQTL data were not used in the MR analysis. All SNPs were identified within 100 kb upstream and downstream of the respective gene locations, with a linkage disequilibrium (LD) threshold of R² < 0.2. The MR analysis was performed using the “TwoSampleMR” package. The MR estimate for each SNP was calculated using the Wald ratio method. If multiple SNPs were available, the MR estimates were weighted by the inverse variance of the ratio estimates (Inverse Variance Weighted, IVW). Cochran’s Q statistic was used to test for heterogeneity, and MR-Egger regression was employed to assess potential horizontal pleiotropy.

Consensus Clustering

After analysis using the “ConsensusClusterPlus” package, unsupervised clustering analysis was performed based on the expression of FRGs, dividing the post-reperfusion samples into different groups. We used t-Distributed Stochastic Neighbor Embedding (t-SNE) to assess sample clusters.

Statistical Analysis

Statistical analyses were performed using R software (version 4.3.0) and GraphPad Prism 9.0. The normality of data distribution was assessed using the Shapiro–Wilk test, and homogeneity of variances was evaluated using Levene’s test. For data that followed a normal distribution and had homogeneous variances, one-way analysis of variance (ANOVA) was used to compare differences among the groups, followed by Tukey’s post-hoc test for pairwise comparisons. If the data did not meet the assumptions of normality or homogeneity of variances, the Kruskal–Wallis H-test was applied, and significant differences were further analyzed using Dunn’s test with Bonferroni correction for pairwise comparisons. Results were presented as mean ± standard deviation of (SD) for normally distributed data or as median and interquartile range (IQR) for non-normally distributed data. Statistical significance was defined as P < 0.05.

Results

FRGs Before and After Reperfusion in Transplanted Pulmonary

Through analysis of the GSE145989 dataset, 483 DEGs were identified, consisting of 423 upregulated and 60 downregulated genes (Figure 1a). Intersection of DEGs with ferroptosis-related genes yielded 24 shared genes, termed FRGs (Figure 1b). The heatmap illustrated the expression of these FRGs in samples before and after lung transplantation reperfusion, revealing upregulation of all FRGs after reperfusion (Figure 1c). These FRGs were categorized into four classes: drivers, markers, inhibitors, and unclassified, comprising 11 drivers, 1 marker, 9 inhibitors, and 7 unclassified (Figure 1d). Some FRGs belonged to multiple classes.

Figure 1 Selection of FRGs and their expressions between samples before and after reperfusion in the GSE145989 dataset. (a) Volcano plot of DEGs. (b) Venn diagram of the intersection between DEGs and ferroptosis-related genes to obtain FRGs. (c) Heatmap of FRGs expressed in samples. (d) Classification of FRGs. Abbreviations: DEGs, differentially expressed genes; FRGs, ferroptosis-related differentially expressed genes.

Functional Enrichment Analysis of FRGs

To investigate the biological functions of FRGs, we conducted GO and KEGG analyses. In the GO enrichment analysis, the top ten enriched biological processes (BP) included cellular response to chemical stress, positive regulation of cytokine production, response to lipopolysaccharide and so on (Figure 2a). Regarding cellular component (CC) and molecular function (MF), FRGs were predominantly associated with cytokine activity, cytokine receptor binding, and molecular sequestering activity (Figure 2a). Furthermore, these FRGs were enriched in KEGG pathways such as the IL-17 signaling pathway and TNF signaling pathway (Figure 2b).

Figure 2 FRGs enrichment analysis results. (a) GO enrichment analysis. (b) KEGG enrichment analyses. Abbreviations: GO, gene ontology; KEGG, Kyoto encyclopedia of genes and genomes; FRGs, ferroptosis-related differentially expressed genes.

Wgcna

In the GSE145989 dataset, WGCNA identified a soft-thresholding power of 4 as optimal (Figure 3a). A total of 9 modules were identified through analysis, among which the turquoise module consisted of 6671 genes and exhibited a strong correlation with lung transplant reperfusion injury (Figure 3b). Additionally, the scatter plot demonstrated a positive correlation (cor=0.813) between the module membership of the turquoise module and gene significance for post-reperfusion lung injury (Figure 3c). Therefore, the turquoise module was considered the key module for subsequent analysis.

Figure 3 Screening shared genes between FRGs and the key module. (a) Selecting soft thresholding power. (b) Heatmap shows the correlation between 9 modules and post- reperfusion. (c) The scatter plot illustrates the correlation between the module membership of the turquoise module and gene significance for post-reperfusion lung injury. (d) The Venn diagram displays the shared genes between FRGs and the turquoise module. Abbreviations: FRGs, ferroptosis-related differentially expressed genes.

Subsequently, the Venn diagram depicted 19 shared genes between FRGs and the key module (Figure 3d).

Identification of Key Genes Using LASSO Regression and SVM-RFE

Additionally, the 19 genes underwent further analysis using LASSO logistic regression (Figure 4a and b) and SVM-RFE algorithm (Figure 4c) to identify key genes. LASSO revealed 9 identified genes (TNFAIP3, KRAS, CXCL2, ZFP36L2, NEDD4L, SESN2, DDIT3, SLC7A5, EZH2), while SVM-RFE identified 11 genes (TNFAIP3, KRAS, CXCL2, GCH1, SLC2A3, MT1G, KDM6B, CDKN1A, SESN2, NEDD4L, IL6). Ultimately, TNFAIP3, KRAS, CXCL2, NEDD4L, and SESN2 were selected as key genes associated with reperfusion (Figure 4d).

Figure 4 Identification of key genes by LASSO regression and SVM-RFE. (a) The tuning parameter (λ) selection plot of LASSO regression. (b) LASSO coefficient profiles of genes associated with post-lung transplant reperfusion. (c) Feature selection with SVM-RFE. (d) Key gene identification: venn diagram. Abbreviations: LASSO, least absolute shrinkage and selection operator; SVM-RFE, support vector machine-recursive feature elimination.

Expression Validation and ROC Curve of Key Genes in the Validation Dataset

To further ascertain the accuracy of key genes as diagnostic biomarkers for lung transplantation ischemia/reperfusion injury, validation was performed using the GSE18995 dataset. The results demonstrated that the expression of all key genes in the validation dataset increased after reperfusion (p<0.05, Figure 5a-e), with ROC curve AUC values all greater than 0.7 (Figure 5f), suggesting that these five key genes may serve as biomarkers for lung transplantation ischemia/reperfusion injury. Furthermore, in the GSE18995 dataset, we conducted a subgroup analysis for donation after cardiac death (DCD) and donation after brain death (DBD) samples. The results revealed that the expression patterns of the key genes were consistent across both the DCD and DBD subgroups, as well as when all samples were combined (data not shown). These findings suggested that the heterogeneity between DCD and DBD samples had a minimal impact on the final conclusions regarding the key genes involved in ferroptosis. This further supported the reliability of these key genes as biomarkers for lung transplantation ischemia/reperfusion injury.

Figure 5 Validation of key genes in GSE18995. (a–e) Changes in expression of key genes before and after reperfusion. (f) ROC curves of key genes. Abbreviation: ROC, receiver operating characteristic.

3.6 Validation of key genes in a mouse model of lung ischemia/reperfusion injury and identification of candidate genes associated with ferroptosis in LIRI.

To determine whether inhibition of ferroptosis during reperfusion could alleviate lung ischemia/reperfusion injury, we administered Fer-1, a ferroptosis inhibitor, intraperitoneally into mice at the onset of reperfusion. HE staining demonstrated that I/R exacerbated lung injury, characterized by neutrophil infiltration, thickening of the alveolar septa, formation of hyaline membranes, and the presence of proteinaceous debris in the alveolar lumen. However, as shown in Figure 6a, Fer-1 ameliorated these morphological abnormalities. Disturbed iron metabolism and lipid peroxidation are critical features of ferroptosis. As depicted in Figure 6b, levels of iron and MDA in mouse lung tissue were significantly elevated following lung I/R, but decreased after Fer-1 treatment. Conversely, GSH, which exerts antioxidant effects, was significantly increased following Fer-1 application. Although the mRNA expression of GPX4 did not show statistical differences among the three groups, there was a trend towards increased expression with Fer-1 treatment. In addition, the mRNA expression of cyclooxygenase 2 (COX2) was elevated after I/R and decreased following Fer-1 application (Figure 6c). These findings suggest that ferroptosis occurs during lung I/R in mice, and that intervention with Fer-1 at the early stage of reperfusion inhibits ferroptosis, thereby attenuating LIRI. The expression of five ferroptosis-associated key genes was further validated in mice with or without lung ischemia/reperfusion injury using RT-qPCR. In mice with lung ischemia/reperfusion injury, the levels of TNFAIP3, CXCL2, NEDD4L, and SESN2 were significantly elevated compared to normal lung tissue, but decreased following Fer-1 treatment. However, the mRNA levels of KRAS did not decrease with Fer-1 application (Figure 6d). Thus, TNFAIP3, CXCL2, NEDD4L, and SESN2 were selected as candidate genes associated with ferroptosis in lung ischemia/reperfusion injury.

Figure 6 Ferroptosis in lung ischemia/reperfusion injury and the levels of ferroptosis-associated key genes. (a) HE staining of lung tissue (×200 magnification). (b) Levels of iron, MDA, and GSH in lung tissue. (c) The mRNA expression of GPX4 and COX2 in lung tissue. (d) The mRNA expression of TNFAIP3, KRAS, CXCL2, NEDD4L and SESN2 in mouse lung tissue. ns: not significant, *: p < 0.05, **: p < 0.01, ***: p < 0.001, and ****: p < 0.0001. Abbreviations: HE, hematoxylin and eosin; MDA, malondialdehyde; GSH, glutathione; GPX4, glutathione peroxidase 4; COX2, cyclooxygenase 2; TNFAIP3, tumor necrosis factor alpha-induced protein 3; KRAS, Kirsten rat sarcoma viral oncogene homolog; CXCL2, C-X-C motif chemokine ligand 2; NEDD4L, neural precursor cell expressed developmentally down-regulated 4-like; SESN2, sestrin 2.

Construction of a Predictive Model and Nomogram for LIRI

Next, we constructed a logistic regression model using four candidate genes (TNFAIP3, CXCL2, NEDD4L and SESN2). The model’s AUC was 0.9891 (95% CI 0.9728–1.0000) (Figure 7a), which was higher than the individual AUCs of the four genes (The AUCs for TNFAIP3, CXCL2, NEDD4L, and SESN2 were 0.9746, 0.9327, 0.8937, and 0.8888, respectively, Figure S1), indicating that the model has better diagnostic ability for LIRI. We also validated the model in the GSE18995 dataset, where the AUC was 0.95 (Figure 7a), which was also higher than the individual AUCs of the four genes, further confirming the robustness of this model. Additionally, we used a nomogram to visualize the risk of disease occurrence. Each variable in the nomogram was projected upward to a point, and the total score of the four variables was converted into an individual’s disease risk. A higher total score corresponded to a higher risk of developing LIRI (Figure 7b). The calibration curve showed no obvious bias between the actual and predicted observations (Figure 7c). We validated the clinical utility of the nomogram using Decision Curve Analysis (DCA). As shown in Figure 7d, the DCA demonstrated a higher overall net benefit across the threshold probability range of 0–1.

Figure 7 A predictive model and nomogram for LIRI. (a) ROC curves showing the diagnostic value of the four-gene model in the GSE145989 and GSE18995 dataset. Blue represents the GSE145989 dataset, while yellow represents GSE18995 dataset. (b) The nomogram predicts LIRI risk based on TNFAIP3, CXCL2, NEDD4L, and SESN2. (c) Calibration curve for nomogram validation. (d) Decision curve analysis based on the nomogram model. Abbreviations: LIRI, lung ischemia/reperfusion injury; ROC, receiver operating characteristic; TNFAIP3, tumor necrosis factor alpha-induced protein 3; CXCL2, C-X-C motif chemokine ligand 2; NEDD4L, neural precursor cell expressed developmentally down-regulated 4-like; SESN2, sestrin 2.

Immune Cell Infiltration

Using the CIBERSORT algorithm, we determined the proportions of 22 immune cell types in the GSE145989 dataset before and after lung transplantation reperfusion (Figure 8a). Compared to before reperfusion, significant changes in immune cells occurred after reperfusion, including Neutrophils, Dendritic cells activated, Mast cells resting, Mast cells activated, Macrophages M2, NK cells resting, Plasma cells, T cells CD4 memory activated, T cells CD4 naive, T cells gamma delta and Macrophages M1 (Figure 8b). Subsequently, correlation analysis was performed on the immune cells, revealing strong synergy between some immune cells, such as Neutrophils and T cells CD4 naive. While, Neutrophils and Macrophages M2 exhibited strong antagonistic effects (Figure 8c). These results suggested that immune cells might play a significant role in influencing lung ischemia/reperfusion injury.

Figure 8 Analysis of immune cell infiltration in the GSE145989 database. (a) Stacked bar plot shows the relative abundance of 22 infiltrating immune cells. (b) Boxplot displays the differences in immune cell infiltration before and after lung transplantation reperfusion. (c) Heatmap shows the correlation between immune cells. (d) Heatmap illustrates the relationship between the candidate genes and immune cells. *: p < 0.05, **: p < 0.01, ***: p < 0.001, and ****: p < 0.0001.

Correlation Between Candidate Genes and Immune Infiltration

Figure 8d illustrated the significant correlations between four candidate genes and immune cells. TNFAIP3 displayed significant correlations with 11 immune cell types, with the most pronounced negative correlation observed with Macrophages M2 (r = −0.698, P < 0.001). CXCL2 showed positive correlation with Neutrophil (r = 0.377, P < 0.001) and the most notable negative correlation with Macrophages M2 (r = −0.544, P < 0.001). NEDD4L demonstrated the most substantial negative correlation with Mast cells resting (r = −0.503, P < 0.001) and the most pronounced positive correlation with Mast cells activated (r = 0.512, P < 0.001). SESN2 was correlated with eight immune cell types, with the strongest correlation observed with M2 macrophages (r = −0.400, P < 0.001). We found a relationship between the candidate genes and immune infiltration, suggesting that these genes may have a potential impact on immune cell infiltration, but whether they are the cause or the consequence remains unclear.

Mendelian Randomization Analysis Between Candidate Genes and PGD

LIRI is the primary cause of PGD following lung transplantation. We utilized Mendelian randomization to further investigate the causal relationship between candidate genes and PGD. Among the four candidate genes, CXCL2 was excluded due to the absence of suitable single nucleotide polymorphisms (SNPs). Figure 9a illustrated the causal effects of each gene on PGD. The inverse variance weighted (IVW) analysis results showed that SESN2 (OR=0.2070, p=2.973×10^-5) was associated with a reduced risk of PGD, whereas the associations of TNFAIP3 and NEDD4L with PGD were not statistically significant. Cochran’s Q test indicated no heterogeneity in the results for SESN2 and MR-Egger regression analysis suggested that there was no horizontal pleiotropy in the results for SESN2. Given the protective role of SESN2 in PGD, we constructed a transcription factor-miRNA network for SESN2 (Figure 9b). The network identified 4 transcription factors and 21 miRNAs that perhaps influence the expression of SESN2.

Figure 9 Mendelian randomization analysis. (a) Mendelian randomization analysis between candidate genes and PGD. (b) The regulatory network construction of SESN2. Abbreviations: PGD, primary graft dysfunction; SESN2, sestrin 2.

Identification of Two Distinct Ferroptosis Patterns

LIRI is a major cause of PGD after lung transplantation, and ferroptosis is thought to exacerbate the pathological process of LIRI. To better understand the clinical heterogeneity of LIRI and its underlying mechanisms, we performed consensus clustering based on 24 FRGs to investigate the relationship between FRG expression and LIRI diversity. After testing various cluster numbers from 2 to 9, we found that two clusters demonstrated the best stability (Figure 10a). And given that some patients develop PGD after lung transplantation while others do not, we thought that classifying LIRI patients into two subgroups better capture this clinical variability. Thus, LIRI patients were categorized into two distinct patterns, Cluster 1 and Cluster 2. T-SNE showed that FRGs could distinguish between the two ferroptosis patterns (Figure 10b). A heatmap was plotted to illustrate the differential expression of the 24 FRGs between the two clusters (Figure 10c). Among them, KDM6B, GCH1, IL6, PROK2, MAP3K14, SLC2A3, MT1G, PTGS2, CDKN1A, PANX1, HMOX1, IL1B, TNFAIP3, SLC7A5, and SESN2 were more highly expressed in Cluster 2 than in Cluster 1, while the differences in the expression levels of KRAS, CXCL2, ZFP36, NEDD4L, DDIT3, EPT1, IFNG, EZH2, and EGR1 between the two groups were not statistically significant. We calculated the ferroptosis score for each sample, with Cluster 2 exhibiting a higher ferroptosis score than Cluster 1 (Figure 10d). Furthermore, ssGSEA was employed to compare the differences in immune cell infiltration between Cluster 1 and Cluster 2 (Figure 10e). We found that Activated B cells, CD56bright natural killer cells, and Immature dendritic cells were significantly more infiltrated in Cluster 1 than in Cluster 2. Effector memory CD4 T cells, Neutrophils, and Type 1 T helper cells showed no statistically significant differences between the two clusters, while the remaining 19 immune cell types were more infiltrated in Cluster 2 than in Cluster 1. To elucidate the potential mechanisms by which these FRGs contribute to the two clusters, we performed GO and KEGG enrichment analyses on the FRGs between the two clusters (Figure 10f). The analysis revealed that functions related to inflammation and immune regulation were significantly enriched. Among these, “regulation of cytokine production involved in immune response” was prominently enriched in the biological process category and intersected with key pathways such as the “TNF signaling pathway” and “IL-17 signaling pathway”, suggesting that the studied targets mediated inflammatory responses through these pathways. Additionally, oxidative stress-related functions were validated across multiple levels. The biological process “response to oxidative stress” was significantly enriched and was closely associated with the molecular function “oxidoreductase activity” and the pathway analysis of the “NF-kappa B signaling pathway”. This indicated that oxidative stress influenced cellular homeostasis through specific signaling pathways. Furthermore, the importance of membrane structures in signal transduction was highlighted. The cellular component terms “plasma membrane raft” and “membrane microdomain” were linked to the molecular function “growth factor receptor binding”, suggesting that these specialized membrane structures served as platforms for signal transmission. These findings indicated that targeting these signaling pathways might alter the subtypes of LIRI. Further studies were required to validate these mechanisms.

Figure 10 Two distinct ferroptosis patterns based on FRGs. (a) Consensus matrices of the 24 FRGs for k = 2. (b) t-SNE plot of the expression profiles of the 24 FRGs. (c) Heat map of the FRGs in Cluster 1 and Cluster 2. (d) Differences in ferroptosis score between Cluster 1 and Cluster 2. (e) Immune cell infiltration in Cluster 1 and 2. (f) GO and KEGG analysis of FRGs in ferroptosis subtypes. ns: not significant, ***: p < 0.001. Abbreviations: FRGs, ferroptosis-related differentially expressed genes; t-SNE, t-distributed Stochastic Neighbor Embedding; GO, gene ontology; KEGG, Kyoto encyclopedia of genes and genomes.

Discussion

In this study, through screening transcriptome data of transplanted lung patients from the GEO database and validating with mouse lung ischemia/reperfusion injury model, we identified and validated several candidate genes related to ferroptosis for LIRI in lung transplantation, including TNFAIP3, CXCL2, NEDD4L, and SESN2. We developed a predictive model based on these genes, demonstrating excellent diagnostic performance. Additionally, through MR analysis, SESN2 was identified as potentially protective against PGD. Establishing a transcription factor-miRNA regulatory network related to SESN2 helped explore the important mechanisms of LIRI. Finally, we identified two distinct subtypes of LIRI by consensus clustering analysis, revealing significant differences related to immune cell infiltration and ferroptosis activity. These findings might serve as diagnostic biomarkers and therapeutic targets for ischemia/reperfusion injury during lung transplantation.

TNFAIP3, also known as zinc finger protein A20, exhibits both deubiquitinating and E3 ubiquitin ligase activities, influencing various signaling pathways including NF-κB and mitogen-activated protein kinase (MAPK), thereby inhibiting the inflammatory pathway and reducing inflammatory responses.^17,18 However, beyond its involvement in inflammation, TNFAIP3 also exerts complex and sometimes contradictory effects on cell death. Recent studies have shown that TNFAIP3 can modulate different types of cell death largely due to its ubiquitin-binding properties, such as limiting necroptosis in T cells and macrophages, promoting autophagy to maintain CD4⁺ T cell survival, and potentially restricting cell survival by promoting apoptosis.^19,20 Regarding ferroptosis, overexpression of TNFAIP3 increased acyl-CoA synthetase long-chain family member 4 (ACSL4) levels, then promoting erastin-induced ferroptosis in human umbilical vein endothelial cells.²¹ In addition, TNFAIP3 in A549 cells treated with erastin increased reactive oxygen species (ROS) production and interacted with ACSL4 and solute carrier family 7 member 11 (SLC7A11) to further promote ferroptosis.²² CXCL2 is a chemokine that can be secreted by monocytes, macrophages, and endothelial cells. It binds to its receptor, C-X-C chemokine receptor type 2 (CXCR2), attracting neutrophils to sites of injury, thereby playing a critical role in acute inflammation and immune responses.²³ CXCL2 has been identified as a key gene in ferroptosis related to ischemia/reperfusion injury in myocardial, renal, and intestinal tissues.^24–26 Consistent with previous findings, our study also found that CXCL2 was positively correlated with neutrophil infiltration in lung tissue after LIRI. NEDD4L, an E3 ubiquitin ligase, plays a vital role in the ubiquitin-proteasome pathway and has been shown to regulate ferroptosis. NEDD4L induced the ubiquitination and degradation of ACSL4, transferrin receptor, CD71, and lactotransferrin, leading to inhibiting ferroptosis by preventing excessive lipid metabolism and iron accumulation.^27–30 Interestingly, recent studies found that NEDD4L could also degrade SLC7A11 and GPX4, promoting ferroptosis.^31–33 Therefore, the role of NEDD4L in ferroptosis is complex, and whether it promotes or inhibits ferroptosis perhaps depends on its specific molecular mechanisms and downstream targets. Our results showed that administering Fer-1 before reperfusion, compared to the LIRI group, alleviated lung injury, decreased levels of iron, MDA, and COX2, and increased levels of GSH and GPX4. This demonstrated that ferroptosis occurs during reperfusion and can be inhibited by Fer-1. Following lung ischemia/reperfusion, the expression of TNFAIP3, CXCL2, NEDD4L were significantly upregulated; however, Fer-1 administration reduced their levels. Additionally, the area under the ROC curve (AUC) of the LIRI prediction model based on these genes far exceeded 0.7, highlighting the potential diagnostic value of TNFAIP3, CXCL2, and NEDD4L in the context of LIRI.

Among the four candidate genes, MR analysis further confirmed that SESN2 was a protective factor against PGD following lung ischemia/reperfusion. SESN2 is a highly conserved oxidative stress protein that protects cells from ferroptosis, with its expression being regulated by the transcription factor Nrf2.^34,35 Conversely, SESN2 could activate Nrf2, which reduced the upregulation of ACSL4 and prevented the downregulation of GPX4, ferritin heavy chain 1, and SLC7A11 in Caco-2 cells under hypoxia-reoxygenation conditions.³⁶ Additionally, SESN2 inhibited ferroptosis in dendritic cells during sepsis by downregulating the activating transcription factor 4 (ATF4)- C/EBP homologous protein (CHOP)- cation transport regulator homolog 1 (CHAC1) signaling pathway and increasing glutathione levels, suggesting SESN2 perhaps protected immune cells from ferroptosis.³⁷ Erastin can significantly increase the expression of SESN2 in a dose- and time-dependent manner.^34,35,38 Our study results showed that SESN2 was elevated after LIRI, and its expression was positively correlated with neutrophil infiltration, possibly as a feedback protective mechanism triggered by LIRI. As a protective factor for PGD, SESN2 can inhibit ferroptosis. Since lung ischemia/reperfusion is a major cause of PGD, we hypothesize that SESN2 reduces the occurrence of PGD by inhibiting ferroptosis caused by ischemia/reperfusion injury during lung transplantation. This hypothesis, however, requires further confirmation through experimental and clinical studies. Additionally, we conducted predictions of transcription factors and miRNAs for SESN2, aiming to enhance its expression and thereby reduce the incidence of LIRI and PGD.

In our study, we observed differential infiltration of immune cells before and after reperfusion, primarily including mast cells, neutrophils, macrophages, dendritic cells, NK cells, and plasma cells. These cells were associated with candidate genes, implying that immune cells might influence the process of ferroptosis, and conversely, ferroptosis might affect immune cell infiltration during lung ischemia/reperfusion. Studies found that extracellular neutrophil extracellular traps (NETs) produced by neutrophils induced ferroptosis in alveolar epithelial cells by promoting m6A methylation of GPX4 and hypoxia-inducible factor 1-alpha (HIF-1α), and also led to ferroptosis of pulmonary endothelial cells.^39–41 Furthermore, macrophages phagocytized damaged or excessive red blood cells, leading to intracellular iron overload, which disrupted the regulatory system of hepcidin and iron transporters, thereby creating conditions conducive to ferroptosis. IL-6, TNF-α, IL-1β produced by macrophages could promoted ferroptosis, while iNOS inhibited ferroptosis.⁴² Additionally, IFN-γ repressed system xc− activation via the JAK/STAT signaling pathway. NK cells can produce IFN-γ, potentially increasing cellular sensitivity to ferroptosis through this mechanism.⁴³ Immune cells might also modulate the genes identified in our study and influence ferroptosis of cells in LIRI through this mechanism. Further studies are needed in the future to validate it.

Personalized medicine is a medical approach that aims to classify diseases accurately and provide personalized treatment based on each individual’s unique characteristics, thereby improving treatment outcomes. In this study, we performed unsupervised consensus clustering analysis on 24 FRGs genes, classifying 67 LIRI patients into two subtypes with different ferroptosis patterns. The results showed significant differences in immune infiltration between these two subtypes, and the ferroptosis score of Cluster 2 was significantly higher than that of Cluster 1. Further GO and KEGG pathway enrichment analysis of the differentially expressed genes indicated that these genes were primarily enriched in pathways related to inflammation and immune regulation, oxidative stress functions, and membrane structures involved in signal transduction. Interventions targeting these signaling pathways may alter the characteristics of LIRI subtypes and provide new therapeutic strategies.

In this study, we screened transcriptome data of transplanted lung patients and validated findings with a mouse lung ischemia/reperfusion injury model, identifying and validating ferroptosis-related genes associated with lung transplantation, including TNFAIP3, CXCL2, NEDD4L, and SESN2. Based on these genes, we developed a predictive model demonstrating excellent diagnostic performance. Additionally, MR analysis identified SESN2 as a potential protective factor against PGD. We also established a transcription factor-miRNA regulatory network related to SESN2, exploring the important mechanisms of LIRI. Finally, through consensus clustering analysis, we identified two distinct LIRI subtypes, revealing significant differences related to immune cell infiltration and ferroptosis activity. These findings may serve as diagnostic biomarkers and therapeutic targets for ischemia/reperfusion injury during lung transplantation.

Despite these significant findings, there are some limitations to our study. Firstly, while we validated the expression of certain genes in a mouse model, their specific roles and relevance in humans require further investigation. Additionally, the hypothesis of SESN2 as a protective factor needs to be confirmed through more extensive experimental and clinical studies to ensure its feasibility and efficacy in clinical applications. By addressing these limitations, we aim to further advance medicine in lung transplantation, providing more effective treatment strategies for LIRI patients.