Identification of Neutrophil Extracellular Trap-Related Gene Expression Signatures in Ischemia Reperfusion Injury During Lung Transplantation: A Transcriptome Analysis and Clinical Validation

Introduction

Lung transplantation has been established worldwide as the only effective remedy for end-stage lung diseases, including chronic obstructive pulmonary disease, interstitial lung disease, cystic fibrosis.^1,2 Over the past few decades, the field of lung transplantation has witnessed remarkable advances, experiencing a substantial surge in lung transplantation procedures since the first successful operation in 1983.³ Annually, over 4000 lung transplantations are conducted worldwide.⁴ Despite, morbidity and mortality of lung transplant patients remain high compared with other solid organ transplants. IRI is an unavoidable occurrence following lung transplantation, which characterized by uncontrolled sterile inflammation, increased microvascular permeability, elevated pulmonary vascular resistance, pulmonary edema, impaired oxygenation, and pulmonary hypertension.^5–7 And it is the main mechanism of primary graft dysfunction (PGD), one of the major clinical challenges for clinicians in the perioperative period of lung transplantation as well as a leading cause of early allograft failure.^8,9 Since the exact underlying pathogenesis of IRI during lung transplantation remains unclear and the absence of effective therapies, further investigation is warranted to explore targeted treatments aimed at the pathogenetically altered expressed genes associated with IRI and PGD.

Neutrophils are the end effectors of innate immune response in IRI after lung transplant.^10,11 Recipient derived neutrophils are recruited to lung allografts and subsequently infiltrate into alveolar space after reperfusion, triggering various damage associated effector events, including the expanded production of ROS, the release of azurophilic granules and metalloproteinase.¹² Among the effector functions of neutrophils, the formation of neutrophil extracellular traps (NETs) generated by a regulated cell death program termed “NETosis” has been confirmed to be closely correlated to PGD following lung transplantation. NETs are composed of extracellular chromatin decorated with histones, antibacterial peptides, and serine proteases.¹³ And those components could contribute to tissue injury and vascular occlusion.^14,15 NETs seemingly function as a biomarker of severity of inflammation. Perioperative elevation of NETs biomarkers has been demonstrated to be associated with PGD3 following lung transplantation, suggesting their potential utility in identifying recipients at risk of PGD3 and initiating preventive therapies.¹⁶ Multiple experimental models of lung transplantation have similarly demonstrated platelet-driven NETs formation in the pathogenesis of PGD.¹⁷ Consequently, the identification of NET-related gene biomarkers will enhance our mechanistic comprehension and offer potential targeted therapeutic strategies for PGD.

Previous transcriptome research of human lung allograft before and after transplantation was largely restricted as the deficiency of clinical data. Therefore, in our study, we used neutrophil infiltration levels to approximately evaluate the severity of IRI after lung transplant based on CIBERSORT algorithm. Subsequently, WGCNA and machine learning algorithms including least absolute shrinkage and selection operator (LASSO) regression and random forest (RF) screened out the most significant genes related to neutrophil infiltration and NETs formation, which may have a potential correlation with PGD. Finally, the expression of hub genes was validated in clinical samples. This is the first time that WGCNA and machine learning were used to explore NETs-related hub genes of IRI after lung transplantation, which could provide novel insight for the early diagnosis at the molecular level and treatment of PGD patients.

Materials and Methods

Data Acquisition and Processing

The flowchart in Figure 1 depicts the methodology of this study. The GSE145989,¹⁸ GSE127003¹⁹ and GSE18995²⁰ microarray data were downloaded from the GEO database (http://www.ncbi.nlm.nih.gov/geo/).²¹ Table 1 presents detailed dataset information. The GSE145989 dataset comprises 67 reperfused lung allograft biopsy samples collected at 1 or 2 hours after reperfusion, extracted on the GPL570 platform (Affymetrix Human Genome U133 Plus 2.0 Array).¹⁸ Our study selectively included 51 reperfused lung allograft biopsy samples collected at 2 hours after reperfusion. All microarray data were normalized through R package “limma”.²² When multiple probes identified the same gene, the average value was calculated to determine its expression. The dataset GSE145989 and GSE127003 were then merged into a larger dataset via R for subsequent analysis. The R package “SVA”²³ was applied to eliminate batch effects by Combat, which was employed to normalize the expression values from different batches or platforms. The GSE18995 dataset was used for validation of hub genes.

Table 1 Basic Information of GEO Datasets Used in the Study

Figure 1 The workflow of the study.

Identification of DEGs and Functional Enrichment Analyses

The normalized and batch-eliminated merged microarray expression dataset underwent differential expression analysis via the R package “limma” to identify differentially expressed genes (DEGs).²² False discovery rate-adjusted p values (FDR) < 0.05 and |log (fold-change)| values ≥0.5 were set as the threshold. The identified DEGs were subsequently subjected to Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analyses. The R package “org.Hs.eg.db.” was initially utilized to convert the hub gene symbols as “ENTREZID”. Subsequently, biological functions were assigned to these genes, including biological process (BP), cellular components (CC), molecular function (MF), and KEGG pathway enrichment, based on the R packages “clusterProfiler”,²⁴ “org.Hs.eg.db”, “enrichplot” and “ggplot2”. A threshold of FDR <0.05 was utilized as the cutoff criterion.

Immune Cell Infiltration Analysis

The relative abundance of 22 immune cell subsets in normalized merged dataset were calculated based on the CIBERSORT gene expression deconvolution package.²⁵ The Wilcoxon test was employed to assess the differences in 22 immune cells between ischemic and reperfused lung allograft biopsy samples. The comparison regarding the proportion of diverse types of immune cells between CIT and reperfusion groups was visualized via the boxplot. A heatmap depicting the correlation of 22 types of infiltrating immune cells was carried out using the R package “corrplot”.²⁶

WGCNA

Since small variation of gene expression data often represents noise, we used Coefficient of Variation values >0.1 (CV genes) to select variant genes in 113 CIT lung samples and 97 reperfused lung samples. Then, the R package “WGCNA” was applied to construct the mRNA co-expression network based on the selected CV genes. Pearson’s correlation matrices were determined with the expression of individual transcripts converted to the similarity matrix. According to the principle of a scale-free network, soft thresholds (power = 5, R² = 0.9) were selected to construct a scale-free co-expression network in turn, and the adjacency matrix was converted into a topological overlap matrix. Then, cluster analysis was performed to identify gene modules, with a minimum size of 30 for the genes dendrogram, and a threshold of 0.25 was set in order to merge identical modules. The module eigengenes were used for component analysis of each module. Pearson’s tests were subsequently used to compute the correlations between the infiltration levels of immune cells and the modules to identify significant modules, with a distinct module with p < 0.05 considered to be significantly correlated with immune cells infiltration. The module with the maximum correlation coefficients was identified as the hub module. Genes in the hub module were further selected with the threshold of membership values (MM) >0.8 and gene significance values (GS) >0.5.

Functional Enrichment Analyses and Construction of PPI Network of Hub Genes

DEGs and brown module genes were intersected to obtain differentially expressed brown module genes through an online Venn diagram tool (http://bioinformatics.psb.ugent.be/webtools/Venn/). Then these genes were used to perform GO and KEGG functional enrichment analyses. The STRING database (https://string-db.org/) was utilized to perform PPI network analysis and identify the interactions among proteins that were encoded by the hub genes with the threshold of interaction score >0.4.²⁷ Then, the results downloaded from the STRING database were visualized through Cytoscape software (v3.9.1).

Machine Learning

Two machine learning algorithms, LASSO and RF, were employed to further filter candidate genes for biomarkers of IRI after lung transplantation based on differentially expressed brown module genes. LASSO, a regression method, serves for variable selection and regularization within statistical models, enhancing both predictive accuracy and comprehensibility.²⁸ The binomial distribution variables are used in the LASSO classification, coupled with lambda. min value for the minimum criterion used to build the model. RF is a suitable approach that offers several advantages, including no restrictions on variable conditions and better accuracy, sensitivity, and specificity, which can be utilized for predicting continuous variables and generating forecasts without apparent variations.²⁹ R packages “glmnet”³⁰ and “randomForest”³¹ were employed for conducting LASSO regression and RF analysis. The intersection genes of LASSO, RF and NET-related genes from literature³² were then considered as hub genes for biomarkers of IRI after lung transplantation.

Differential Expression Analysis and ROC Curve Validation

R packages “limma” and “ggpubr” were used to analyze the differential expression levels of the final hub genes in CIT and reperfusion groups based on merged dataset and GSE18995. Meanwhile, receiver operating characteristic (ROC) curve analysis was performed for each hub gene using the R package “pROC”,³³ and the area under the curve (AUC) with a 95% confidence interval (CI) was calculated. The significance of IRI was judged by the AUC value, with values close to 1 indicating higher accuracy of the model training.

Construction of Nomogram and Verification

Based on the four candidate genes, the R package “rms” was applied to construct the nomogram. “Points” indicates the score of candidate genes, and “Total Points” signifies the summation of all the scores of genes above. The ROC was subsequently established to assess the predictive value of the nomogram and candidate genes regarding IRI after lung transplantation. Calculation of AUC and 95% CI quantified their predictive capacity. AUC > 0.7 was considered as ideal predictive value. Then calibration curve was plotted to assess the predictive accuracy of the nomogram. Furthermore, the decision analysis (DCA) and clinical impact curve were plotted to evaluate the clinical value of the nomogram.

Construction of Potential TF-Target Gene Regulatory Networks and Small-Molecule Drug Prediction

DAVID database (https://david.ncifcrf.gov/) was employed to predict upstream transcription factors (TFs) based on DEGs. The results were then visualized using the Cytoscape software. The Broad Institutes Connectivity Map (CMap) database (https://clue.io/) was employed to identify potential small candidate molecules. We queried CMap with a signature containing the top 150 most upregulated genes in DEGs and obtained an inversely correlated connectivity score (CS) to IRI. Four small molecules selected according to CS could reverse biological states encoded in specific gene expression markers, which represents that these compounds might have a potential therapeutic effect in IRI after lung transplantation. Moreover, we obtained detailed information, including targeted genes, mechanism of action (MOA), and clinical phase, as well as 3D confirmation of the established small molecules extracted from PubChem (https://pubchem.ncbi.nml.nih.gov).

Clinical Samples Collection and Treatment

A total of 12 patients who underwent lung transplantation with PGD 2 grade or more, as well as 12 patients with no PGD, were selected and enrolled into this study from Shanghai Pulmonary Hospital between January 2022 and August 2022 according to 2016 consensus statement of ISHLT working group on PGD definition.³⁴ Patients were stratified by the cut point of PaO2/FiO2 ratio of 300 at 48 hours after reperfusion and the presence of bilateral alveolar infiltrates on chest X ray. The use of postoperative extracorporeal life support (ECLS) should be explicitly reported and taken into account. Furthermore, utilizing atomized prostacyclin or other drugs that may improve oxygenation did not affect PGD classification.³⁵ Human peripheral blood samples at cold ischemic time and 12 hours after reperfusion were obtained from these lung transplant patients, collected into citrated vacutainers and centrifuged at 2000 ×g for 20 min. And then the generated plasma was separated and preserved at −80°C for further experiments. The study was approved (No. K22-267Z) by the Research Ethics Commission of Shanghai Pulmonary Hospital (Shanghai, China), and written informed consent was obtained from the recipients involved in the study or their representatives.

Neutrophils Isolation and Immunofluorescence Staining

Human peripheral neutrophils were isolated by discontinuous density gradient centrifugation with two solutions of differential density (Ficoll and Histopaque-1119, Sigma; St Louis, MO, USA), according to the manufacturer’s instructions. NET formation was examined by immunofluorescence staining. The Human peripheral neutrophils were first treated with SYTOX Green (Thermo Fisher Scientific, Waltham, MA, USA) for 15 min. Subsequently, the neutrophils were fixed in 4% paraformaldehyde and permeabilized with Triton X-100. Then, fixed neutrophils were stained with rabbit polyclonal anti-citrullinated histone H3 (CitH3) rabbit antibody (1:1000 dilution; Abcam, Cambridge, UK) overnight at 4°C. These slides were incubated with goat anti-rabbit IgG labeled-secondary antibodies (EarthOx, San Francisco, CA, USA) for 1 h and the DNA was stained with 4ʹ,6-diamidino-2-phenylindole (DAPI) at room temperature for 10 min. Last, the cells were washed twice with PBS and covered with an antifluorescence quencher. A fluorescence microscope was used to image NETs to assess their structure and location.

Quantification of Plasma dsDNA

We used the PicoGreen dsDNA Quantification Kit (Solarbio, Beijing, P.R. China) to quantify plasma NETs levels. All of the experimental manipulations were based on the instructions of manufacturers with minor adjustments. Each sample was analyzed twice.

Enzyme-Linked Immunosorbent Assay

Plasma levels of three proteins (S100A12, PADI4, and MMP9) were quantified using enzyme-linked immunosorbent assay (Cusabio Biotech Co., Wuhan, Hubei Province, P.R. China). ELISA kits for S100A12 (CSB-E13095h), PADI4(CSB-E16219h) and MMP9 (CSB-E08006h) were purchased from Cusabio, Wuhan, China (https://www.cusabio.com/). All of the experimental manipulations were based on the instructions of manufacturers with minor adjustments.

Statistical Analysis

All statistical analyses were conducted using R (version 4.2.1) or GraphPad Prism 9.5.1. Intergroup comparisons were performed by Student’s t-test for normal distributed variables and Mann–Whitney test for skewed distributed data. Correlations were assessed through Spearman correlation or Pearson correlation analysis. *p<0.05, **p<0.01, and ***p<0.001 were considered statistically significant.

Results

Data Preprocessing and DEG Screening

The dataset GSE145989 and GSE127003 downloaded from GEO were merged and batch-normalized for subsequent analysis, resulting in a larger merged dataset with 113 ischemic lung allograft biopsy samples and 97 reperfused lung allograft biopsy samples. A total of 586 DEGs were identified using the merged dataset based on the criteria of |logFC| > 0.5 and adjusted p< 0.05, with 457 genes up-regulated and 129 genes down-regulated presented in the volcano plots (Figure 2A). The top 30 genes from the upregulated and downregulated genes were selected to create a heat map, where red represents upregulated genes and blue represents downregulated genes (Figure 2B). GO and KEGG functional enrichment analysis was performed on the DEGs with top 10 enrichment terminologies presented in Figure 2C and D. The GO analysis revealed predominant enrichment of these genes in cytokine−mediated signaling pathway, leukocyte migration, receptor ligand activity, signaling receptor activator activity, chemokine activity, and chemokine receptor binding (Figure 2C). The KEGG pathway enrichment analysis identified the top four significantly enriched pathways as TNF signaling pathway, cytokine−cytokine receptor interaction, IL−17 signaling pathway, and NF−kappa B signaling pathway (Figure 2D). These results strongly indicated the involvement of immune response, especially in leukocyte migration and inflammatory response, in the development and occurrence of IRI during lung transplantation.

Figure 2 DEGs enrichment analysis. (A) Volcano plot of 586 DEGs. The X-axis represents log2 FC, and the Y-axis represents the log-transformed adjusted P values. (B) Heatmap of the top 30 DEGs. (C) GO analysis of DEGs: biological process (BP), cellular composition (CC), and molecular function (MF) of hub genes. (D) KEGG pathway analysis of DEGs revealed the top 10 relevant pathways: The dot sizes refer to the number of genes involved in the pathway, and the dot colors depict the p-values.

Estimation of the Immune-Infiltration Level in IRI After Lung Transplantation and Construction of Weighted Gene Co-Expression Network

CIBERSORT, a bioinformatics algorithm capable of estimating 22 immune cell types, was utilized to assess the composition of immune cells based on the gene expression matrix.

Using the CIBERSORT algorithm, the analysis evaluated the infiltration levels of 22 distinct types of immune cells in the ischemia and reperfused groups, presented in the boxplot (Figure 3A). The results of immune infiltration analysis indicated an increased abundance of immune cells in the allograft after reperfusion, including neutrophils, CD4+ T cells, NK cells, monocytes, dendritic cells, and mast cells. Especially, the abundance of neutrophil in the lung allograft exhibited a significant increase following reperfusion. The correlation of 22 types of immune cells revealed that neutrophils were positively associated with resting NK cells (r = 0.43) and Monocytes (r = 0.27), whereas negatively related to macrophages M1 (r = −0.48) and macrophages M2 (r = −0.58) (Figure 3B). The extent of immune cell infiltration in the allograft can be indicative of inflammatory damage to the lung allograft, especially the neutrophils serve as the end effector of IRI, and this was employed as a trait data for the subsequent construction of WGCNA.

Figure 3 CIBERSORT immune cell infiltration analysis and WGCNA to cluster genes into different modules. (A) Comparison of the scores for immune cells estimated by the CIBERSORT algorithm between the CIT and reperfusion groups. *p<0.05, **p<0.01, ***p<0.001, ****p<0.0001. (B) The correlation analysis of 22 types of immune cells in lung allograft. (C) Estimation of the scale Independence index of the 1–20 soft threshold power (β = 5) and determination of the mean connectivity of the 1–20 soft threshold power. (D) Dendrogram of all CV genes. (E) Relationships of consensus module eigengenes and macrophages. The rows in the figure correspond to consensus modules, and the columns correspond to macrophage subtypes. The numbers in each module indicated the correlation coefficients to show the association between the corresponding module and macrophages, along with the p values shown below in parentheses (red indicates a positive correlation and blue indicates a negative correlation). (F) Module membership in brown modules and gene significance for neutrophil. Each dot indicates a gene. Genes with module membership values > 0.8 and gene significances > 0.5 inside the red box are candidate hub genes.

To further precisely identify the central hub genes associated with IRI after lung transplantation, we constructed a gene co-expression network using the WGCNA algorithm. A total of 3779 CV genes were recognized for the construction of WGCNA (Supplementary Table S1). We choose immune cell infiltration extent in the CIBERSORT results as the trait data for WGCNA construction (Supplementary Table S2). The soft threshold was set to 5 to satisfy the scale-free topology of the network, resulting in an R² of approximately 0.9 and notably high average connectivity (Figure 3C). Using the average linkage hierarchical clustering method, we constructed a gene hierarchy clustering dendrogram, identifying 9 similar gene modules (Figure 3D). Among these 9 modules, the correlation (R²) of the brown module with neutrophils notably exceeded that of the other modules, reaching 0.73 (Figure 3E). Consequently, the brown module was identified as the hub module. The scatter plot (Figure 3F) showed a strong correlation between GS and MM in the brown module (Coy=0.91, p=3.2e-110). A total of 75 hub genes with module membership values >0.8 and gene significance values >0.5 were discerned from the brown module (Supplementary Table S3).

Functional Enrichment and PPI Network Analysis of the Hub Genes

For further exploration of IRI mechanism during lung transplantation, we intersected DEGs and brown module genes and obtained a total of 38 signature genes illustrated in Figure 4A (Supplementary Table S4). Subsequently, these signature genes were applied GO and KEGG pathway enrichment functional analysis (Figure 4B and C). GO analysis showed that these genes were mainly enriched in leukocyte migration, positive regulation of response to external stimulus, leukocyte chemotaxis, and cell chemotaxis, immune receptor activity, and cytokine receptor activity (Figure 4B). The KEGG pathway analysis indicated significant enrichment of the hub genes primarily in apoptosis, neutrophil extracellular trap formation, and cytokine receptor, and cytokine–cytokine receptor interaction pathways (Figure 4C). The protein–protein interaction (PPI) network was constructed to further assess the interactions among these genes using the STRING online tool (https://stringdb.org/) (Figure 4D). The Cytoscape software (v3.9.1) was applied to visualize the PPI network, comprising 28 nodes. All of these 38 hub genes were selected for subsequent analysis.

Figure 4 Identification and enrichment analysis of signature genes. (A) Venn diagram of DEGs and hub genes in brown module. Commonly intersecting genes were identified as signature genes. (B) GO analysis of signature genes. (C) KEGG analysis of signature genes. (D) Protein–protein interaction (PPI) network of signature genes in was constructed using the STRING database.

Identification of NET Formation-Related Biomarkers via Machine Learning

The LASSO regression and RF machine learning algorithms were applied to further streamline the important characteristic variables for the identification of hub genes. The LASSO regression algorithm identified 34 potential candidate biomarkers based on lambda. min (Figure 5A and B). The RF algorithm in combination with feature selection was used to determine the association between error rate, number of classification trees, and the top 25 genes by weight were selected (Figure 5C and D). Four genes (S100A12, DAPI4, FCAR, and MMP9) were ultimately obtained by intersecting the genes screened from the LASSO and RF algorithm, alongside 136 NETs related genes from literature, visualized represented in Venn diagram (Figure 5E). The four key genes (S100A12, DAPI4, FCAR, and MMP9) were identified as the final hub genes. Gene correlations, depicted in Figure 5F, indicated a positive correlation among the four genes, suggesting their substantial functional similarity.

Figure 5 Identification of final four key NETs related genes. (A and B) The candidate hub genes obtained by LASSO regression with 10-fold cross-validation. 34 potential candidate biomarkers were identified based on lambda. min. (C and D) The randomForest error rate versus the number of classification trees; and 28 genes are ranked based on the importance score. (E) Venn diagram shows that four final genes are identified via the above two algorithms as well as 136 NETs related genes. (F) Correlation between the four final hub genes.

Differential Expression Analysis and Validation of Hub Genes

The four final key genes were performed differential expression analysis based on the merged data (Figure 6A). Compared with CIT groups, the expression of the four final hub genes was upregulated in reperfusion groups (p < 0.001) (Figure 6A). Subsequently, further validation of the four final hub genes was conducted using the external dataset GSE18995. Similarly, all four hub genes exhibited higher expression levels in the reperfusion groups compared to the CIT groups (Figure 6B). The ROC curve analysis was utilized to evaluate the sensitivity and specificity of the four hub genes in external dataset for IRI assessment, revealing AUC values exceeding 0.70 for all four genes, indicating their significant representation for IRI (Figure 6C).

Figure 6 Validation of the diagnostic values of the key genes. (A) Differential expression of four key genes between the ischemic and reperfused lung allografts based on the merged dataset. ***p<0.001. (B) Differential expression of four key genes between the ischemic and reperfused lung allografts based on external dataset GSE18995. ***p<0.001. (C) ROC curves for the validation of FCAR, MMP9, PADI4, and S100A12 based on external dataset GSE18995.

Construction and Assessment of a Nomogram for IRI

The R package “rms” was applied to establish a nomogram for IRI after lung transplantation based on the four candidate hub genes (Figure 7A). ROC curves were established to assess the predictive specificity and sensitivity of both the nomogram model (Figure 7B) and each individual gene (Figure 7C). The ROC curve analysis showed that the AUC value of the nomogram model was higher than any other candidate gene, signifying the accuracy of the nomogram model. Compared with single gene, this nomogram achieved a better predictive value. Additionally, the calibration curve confirmed the accuracy of the nomogram model by demonstrating minimal error between the actual and predicted risk (Figure 7D). DCA manifested that the “nomogram” curve exceeded the gray line, and the “PADI4, MMP9, S100A12, and FCAR” curve implied that patients could benefit from the nomogram model at a high-risk threshold of 0 to 0.9. The nomogram model provided a greater clinical benefit than the “PADI4, MMP9, S100A12, and FCAR” curve (Figure 7E). The clinical impact curve was plotted to evaluate the clinical value of the nomogram model (Figure 7F). These results implied that the four genes may play a key role in the process of IRI after lung transplantation.

Figure 7 Construction and validation of a nomogram model for IRI diagnosis. (A) The nomogram of diagnostic biomarkers to predict the occurrence of IRI. (B) The ROC curve of nomogram model to predict IRI diagnostic value. (C) The ROC curve of each candidate hub gene (PADI4, MMP9, S100A12, and FCAR) to predict IRI diagnostic value. (D) The calibration curve to assess the predictive power of the nomogram model. (E) The DCA curve to evaluate the clinical application value of nomogram model. (F) Clinical impact curves of the nomogram model.

Construction of Potential TF-Target Gene Regulatory Networks and Identification of Small Molecular Therapeutic Agents

The interaction network consisted of 82 genes and 2 TFs (Figure 8A). NFKAPPAB was identified to regulate 50 genes, including S100A12, and PADI4. AP1 could regulate 50 genes, including PADI4 and EREG. The CMap analysis was applied to search therapeutic compounds of which gene expression patterns were opposite to IRI. Figure 8B shows that 4 small molecules with the highest absolute CMap scores were eventually identified as therapeutic agents for IRI after lung transplantation. Their detailed information, including targeted genes, MOA, and clinical phase were shown in Table 2.

Table 2 Detailed Information of the 4 Small Molecules with the Highest Absolute CMap Scores

Figure 8 Construction of TF-Target Gene Regulatory Networks and Identification of small molecular therapeutic agents. (A) TF-target gene regulatory network of DEGs. (B) 3D conformers of the top 4 candidate therapeutic agents with the highest absolute CS scores. The DEGs and Broad Institutes Connectivity Map (cMAP) database (https://portals.broadinstitute.org/cmap) were used for the identification of small molecules. The threshold CS score chosen was set at < −95. Molecules displaying the most elevated absolute CS scores indicate gene expression pattern that stand in opposition to the specific expression patterns observed in IRI. This observation implies their promising therapeutic potential for mitigating IRI subsequent to lung transplantation.

Baseline Characteristics of Lung Transplant Patients

The baseline characteristics of patients enrolled in this study are shown in Table 3. No differences were observed between groups concerning age, sex, body mass index (BMI), transplant diagnosis, pretransplant mean pulmonary artery pressures (mPAP), or transplant type. Patients with a PGD grade of 2 or more exhibited longer durations of mechanical ventilation compared to non-PGD patients.

Table 3 Preoperative Characteristics of Lung Transplant Patients Stratified by Subsequent Development of Grade 2 or Above Primary Graft Dysfunction

Plasma NET Formation Increased in Ischemia-Reperfusion During Lung Transplantation

Neutrophil immunofluorescence staining and plasma dsDNA quantification were employed to confirm the occurrence of NETs formation during ischemia-reperfusion in lung transplantation. The immunofluorescence image showed that NET component H3cit colocalized to extracellular chromatin fibers of neutrophils in reperfusion groups but hardly in CIT groups (Figure 9A). The levels of plasma dsDNA were significantly higher in reperfusion groups (Figure 9B). These results confirmed NETs formation in ischemia-reperfusion during lung transplantation. Furthermore, compared with non-PGD patients, both pre-reperfusion and 12 hours after reperfusion levels of plasma dsDNA in patients with PGD grade 2 or more were significantly higher (Table 4).

Table 4 Mean MMP9, PADI4, dsDNA, and S100A12 Expression Levels in Sera of Lung Transplant Patients Stratified by Development of Grade 2 or Above Primary Graft Dysfunction

Figure 9 Validation of clinical samples. (A) Immunofluorescence microscopy showing neutrophil extracellular traps, defined as SYTOX Green (green), histone (red), and colocalized extracellular DNA (blue), present in lungs allograft after reperfusion. Scale bar = 130 μm. Results are representative of at least three independent experiments. (B) The levels of dsDNA in the plasma of patients undergoing lung transplantation. *p<0.05, ***p<0.001. (C–E) The levels of MMP9, PADI4, and S100A12 in the plasma of patients undergoing lung transplantation. *p<0.05, ***p<0.001. (F) Correlation analysis between dsDNA and MMP9 expression levels. (G) Correlation analysis between dsDNA and PADI4 expression levels. (H) Correlation analysis between dsDNA and S100A12 expression levels.

Validation of the Hub Gene Expression in Clinical Samples via ELISA

To verify the changes in the expression levels of hub genes during the process of lung ischemia-reperfusion, we detected the expression of MMP9, PADI4, and S100A12 in plasma through ELISA. The results showed that MMP9, PADI4, and S100A12 expression was significantly upregulated after reperfusion in total of 24 patients underwent lung transplant (Figure 9C–E). Specifically, among 12 patients with PGD 2 grade or more, a significant increase of PADI4 and S100A12 expression was observed, whereas the MMP9 and S100A12 levels in plasma increased among 12 non-PGD patients after lung transplant. In addition, MMP9 and S100A12 had higher pre-reperfusion levels in patients with PGD 2 grade or more and only 12-hours S100A12 levels are significantly higher in patients with PGD 2 grade or above compared with non-PGD patients. Changes in dsDNA correlated with the changes in MMP9 and S100A12, yet no significant correlation was found between changes in dsDNA and PADI4 (Figure 9F–H).

Discussion

In the current study, we identified four novel NET-related biomarkers during lung ischemia-reperfusion through transcriptome analysis, and further validated their expression in clinical samples. NETs, considered to be web-like structures composed of DNA and granule proteins released after cell death, were found induced in mouse models of pulmonary IRI and human lung allograft after transplantation, with the pathogenicity in PGD.^16,17 NETs released from activated neutrophils contribute to inflammatory process and already formed NETs further recruit neutrophils, creating a feedforward cycle.³⁶ Elevated perioperative NETs biomarkers, including citrullinated H3R8 nucleosome histones and circulating H3.1 nucleosomes, indirectly measuring NET formation, have been demonstrated to be related to PGD grade 3.^16,37 Considering the particular sensitivity to DNAse-1 mediated digestion,³⁸ researchers attempted to attain the breakdown of NETs by DNAse-1 to restore lung function. Despite the rapidly degraded NETs in allograft by DNAse-1 treatment, the consequential release of NET fragments contributes to consistent interactions between infiltrating CD4+ T cells and donor‐derived antigen presenting cells as well as the activation of toll‐like receptor (TLR)-myeloid differentiation factor 88 (MyD88) signaling pathways, which subsequently promotes human alveolar macrophage inflammatory cytokine production.³⁹ A recent study showed the use of a cytokine adsorber during lung transplantation reduced NETs and the incidence of PGD at 1-month post-transplant, indicating the availability of NETs reduction in alleviating PGD.³⁷ Nevertheless, gene-targeted blocking of NET formation presents another highly valuable approach in mitigating lung IRI.

Further analysis via LASSO and RF identified 4 key genes significantly associated with NET formation: FCAR, MMP9, PADI4, and S100A12. All of these key genes upregulated in the lung allograft 2h after reperfusion, indicating the initiation of the programmed NETosis launched in the early phase after reperfusion. Elevated plasma dsDNA levels in lung transplant patients indicated abundant NET formation, with severe NETs formation observed in patients with PGD 2 grade or above. Interestingly, higher pre-reperfusion dsDNA levels were observed in patients with a PGD grade of 2 or higher, indicating a potentially elevated pre-transplant NETosis degree and inflammatory status correlating with PGD incidence.

In addition, we applied CMap database to identify distinct small molecules effective in IRI. Mitoxantrone, approved for treating aggressive relapsing and progressive multiple sclerosis, is a topoisomerase inhibitor exhibiting anti-inflammatory and antioxidant properties.⁴⁰ PD-98059, an inhibitor of the MAPK/ERK signaling pathway, plays a pivotal role in cell survival and inflammation regulation.⁴¹ PD-98059 confers a protection against ERK-mediated apoptosis.⁴² Nevertheless, its involvement in NETosis remains ambiguous. Triptolide, derived from Thunder God Vine, is recognized for its anti-inflammatory, immunosuppressive, and antioxidative properties.⁴³ Vemurafenib is a B-RAF inhibitor primarily used for treating malignant melanoma.⁴⁴ Although its primary application is in cancer therapy, studies indicate the involvement of the B-RAF pathway in inflammation and oxidative stress.⁴⁵ IRI in lung transplantation represents a multifaceted and intricate pathology, characterized by the intricate interplay of oxidative stress, metabolic disruptions, immune activation, and cellular demise.⁴⁶ Hence, our discoveries emphasize the potential benefits of targeting specific pathways or receptors to regulate immune cell function and mitigate inflammation in various aspects of lung transplantation.

FCAR, encoding FcαRI (CD89), the human Fc receptor for IgA, is highly expressed on neutrophils.⁴⁷ Previous studies demonstrated that FCAR shapes neutrophil response and innate immunity, inducing rapid NET formation and NETosis by IgA-opsonized particles in the presence of bacteria.⁴⁸ Additionally, a study demonstrated the upregulation of FcαRI surface expression in neutrophils following exposure to inflammatory mediators like GM-CSF, the TLR4 agonist LPS, and the cytokines G-CSF and TNF-α,⁴⁹ which is consistent with the upregulation of FCAR after reperfusion during lung transplantation in our current study. However, Wehrli et al previously reported that FcαRI induces distinct forms of neutrophil death, depending on the inflammatory microenvironment. FcαRI-mediated neutrophil death is apoptotic in the absence of bacteria, but under GM-CSF or LPS stimulation it still includes the partially nonapoptotic death, involved caspases, ROS, PI3K, and MAPKs.⁴⁹ The observed upregulation of FCAR in our study suggests a potential association with neutrophil activity, including NETs formation. Further research is necessary to elucidate the role of FCAR in neutrophil response and NETs formation during the sterile inflammation of IRI after lung transplantation.

S100 calcium-binding protein A12 (S100A12), a member of the S100 family of EF-hand calcium-binding proteins, is predominantly expressed and secreted by activated neutrophils, accounting for approximately 5% of the total amount of cytosolic proteins in neutrophil granulocytes.⁵⁰ Extracellular S100A12 shows cytokine-like characteristics and serum S100A12 serves as a marker of various inflammatory diseases.^51–53 Various clinical evidence demonstrated that S100A12 can be markedly upregulated in diverse ischemic inflammatory conditions, suggesting its potential as a sensitive and specific diagnostic and prognostic marker of ischemic diseases.^54,55 Hongxing Lei recently found that activation of S100A12 was observed in COVID-19 patients with its activation levels correlated with COVID-19 severity based on transcriptome datasets derived from the blood of COVID-19 patients.⁵⁶ PGD shares similar pathological features as COVID-19, characterized by uncontrolled inflammatory response leading to diffuse alveolar injury and progressive hypoxemia. Consistent with COVID-19 results, our transcriptome results indicated a significant upregulation of S100A12 expression following lung ischemia-reperfusion, as well as the elevated S100A12 plasma levels associated with PGD progression. Mechanistically, S100A12 interacts with the multiligand receptor for advanced glycation end products (RAGE) and toll-like receptor 4 (TLR4), pivotal in activating inflammatory cells like macrophages and lymphocytes, thereby promoting inflammation.^57,58 Animal models showed that S100A12 promotes inflammation and cell apoptosis in sepsis-induced ARDS via activation of NLRP3 inflammasome signaling.⁵⁹ Zhang et al demonstrated that S100A12 promotes inflammation and apoptosis in ischemia/reperfusion injury via ERK signaling through in vitro experiments.⁶⁰ These findings suggest a potential role for S100A12 in lung ischemia-reperfusion, and plasma levels might serve as biomarkers of PGD progression, requiring further mechanistic investigation.

Proteins of the matrix metalloproteinase (MMP) family can degrade almost all extracellular matrix (ECM) components derived from inactive precursor zymogens.⁶¹ MMP9, a representative member of MMPs, degrades the ECM components of the basement membrane to allow neutrophil migration in the early phase of IRI. Yano et al observed elevated MMP9 activity and gene expression in IRI using rat lung transplantation models consistent with our findings in human lung transplantation transcriptome analysis, as well as a decreased level of tissue inhibitors of metalloproteinase (TIMP)-1, which is ubiquitous natural inhibitor and form complexes with MMP-9.⁶² The MMP9 mRNA to TIMP-1 mRNA ratio peaked in the early reperfusion phase, notably higher than in other phases, suggesting its potential as an early biomarker for IRI.⁶³ In addition, increased activity of Caspase-1 and IL-1β were observed after IRI, and the lung pyroptosis and injury alleviated when treatment the mice with MMP9 selective inhibitor SB-3CT, suggesting MMP9 involvement in IRI after lung transplant, potentially associated with promoting lung pyroptosis.⁶⁴ The transcriptome analysis in this study showed elevated MMP9 mRNA levels after reperfusion, and we observed a significant difference in MMP9 plasma levels between pre- and post-transplant of total 24 lung transplant patients. Nonetheless, no significant peri-reperfusion differences were observed in patients with PGD 2 grade or above, possibly due to the limited sample volume. The correlation analysis showed the changes in plasma MMP9 levels were correlated with NET formation, but the specific mechanism in PGD requires further research.

Peptidyl arginine deiminase 4 (PAD4), encoded by PADI4 and predominantly expressed in granulocytes, is a crucial post-translational modifying enzyme that converts arginine residues into citrulline residues in the presence of calcium ions. It has been demonstrated that PAD4-mediated histone citrullination via converting arginine to citrulline is the prerequisite and trigger for NETs formation. This process reduces charge-based interactions with DNA, promoting chromatin decondensation.⁶⁵ The observed positive correlation of plasma PAD4 levels and plasma dsDNA levels in our study future indicates the potential relevance between PAD4 and NETs. NET formation involves two mechanisms: neutrophil elastase and other proteases in granules possibly cleaving histones, and PAD4 converting arginine to citrulline, leading to histone release from DNA, essential for NET release.⁶⁶ Studies have revealed that efficient DNA decondensation, nuclear lamin meshwork, and nuclear envelope rupture and extracellular DNA release in NETosis require PAD4 enzymatic activity and nuclear localization signal.⁶⁷ Moreover, animal studies showed that PADI4 deficiency could ameliorate NET formation and inflammatory response induced by NETs.^68,69 Our transcriptome analysis revealed a significant increase in PADI4 expression, particularly in patients with PGD 2 grade or more, and the upregulation of PADI4 may contribute to the elevated NET formation. Nevertheless, our results indicate no significant alteration in PADI4 expression levels after reperfusion in non-PGD patients. Correlation analysis showed a deficient correlation between PADI4 and dsDNA levels, possibly influenced by the baseline status in lung transplant patients. Further multi-central and high-volume clinical research could resolve this issue. In addition, higher post-reperfusion PADI4 expression in plasma of patients with PGD grade 2 or more compared with non-PGD patients indicates that NET formation is a vital segment of the pathogenesis of IRI in PGD. Therefore, prevention of NET formation by inhibition of PADI4 to block histone release maybe an attractive strategy to prevent NET‐induced tissue damage in lung transplantation.

However, our study still has a few limitations. First, the limited sample size in our validation study might restrict the generalizability of the results. Therefore, larger-sample studies and functional studies of key genes in lung IRI are warranted. Besides, the final hub genes just were validated at the protein level and a PCR validated at the gene level using allograft tissue is required. Moreover, in addition to clinical samples, animal models of IRI and lung transplantation are necessary for the further validation of final hub genes. Future research requires further exploration of the mechanisms underlying these hub genes based on larger size clinical samples and animal models. Lastly, the NETs related genes were based on current literature, requiring continuous exploration and updates. The limitations of our research also require attention.

Our research provides evidence for peripheral blood biomarkers to predict inflammatory injury after lung transplantation. More importantly, peripheral blood samples are easily obtainable, repeatable, and minimally invasive for lung transplant patients. The results of the current study may catalyze further exploration of immune-associated mechanisms in lung transplantation in the future. With the utilization of bulk and single-cell RNA sequencing in broadening research as the decreasing cost, more specific factors which remarkably affect transplant outcomes will be allowed to identify.

Conclusion

In conclusion, we screened out four final hub genes as NETs-related biomarkers of IRI during lung transplantation based on WGCNA and machine learning algorithms. S100A12, PADI4, and MMP9 were validated to be highly expressed after reperfusion via clinical samples, and have been further confirmed to be related with the development of PGD. Our study provides novel insight into IRI during lung transplantation at the immune and molecular levels; nevertheless, the potential NETs related biomarkers identified in our study require more laboratory data for further validation.