Identification and Validation of Hub Genes Related to Neutrophil Extracellular Traps-Mediated Cell Damage During Myocardial Infarction

Introduction

Myocardial infarction (MI), caused by thrombosis or arterial occlusion, is the leading cause of morbidity and mortality among all cardiovascular diseases.¹ Many studies have shown that immune cells are involved in the inflammatory response to MI and heart failure.^2–5 For example, myocardial injury triggers an infiltration of neutrophils and macrophages into the heart tissue. These immune cells play a crucial role in removing cellular debris, initiating inflammation through the secretion of pro-inflammatory cytokines, and attracting other pro-inflammatory cells.^2,6

Neutrophils, the most predominant type of leukocyte, arrive at the site of inflammation in a cascade-like manner. This results in the activation of specific effector functions such as degranulation, phagocytosis, and the release of reactive oxygen species.⁷ In 2004, Brinkmann et al demonstrated that neutrophils stimulation with phorbol myristate acetate (PMA) release granule proteins and chromatin that together form extracellular fibers that bind Gram-positive and -negative bacteria.⁸ Subsequent studies by Fuchs et al showed that upon stimulation, the nucleus of neutrophils loses its shape, followed by disintegration of the nuclear envelope and granule membranes, allowing for the mixing of neutrophil extracellular traps (NETs) components, and finally, NETs are released with the rupture of the cell membrane.⁹ Distinct from apoptosis and necrosis, this cell death process is dependent on the creation of reactive oxygen species (ROS) by NADPH oxidase.⁹ This will eventually lead to the release of myeloperoxidase (MPO) and neutrophil elastase (NE) from cytoplasmic azurophilic granules.¹⁰ Except for degrading virulence factors and killing bacteria, there is growing evidence that NETs play a role in a variety of diseases, such as sepsis,¹¹ acute respiratory distress syndrome (ARDS),¹² chronic obstructive pulmonary disease (COPD),¹³ systemic lupus erythematosus (SLE),¹⁴ and atherosclerosis.¹⁵ Recently, many studies have shown that NETs play a damaging role in acute myocardial infarction (AMI).¹⁶ For example, one study found that in addition to providing a scaffold for thrombus formation, NETs formation also exacerbates endothelial cell injury, which may be an important contributor to the development of MI.¹⁷

Here, we conducted a comprehensive analysis of NETs-related genes in MI, with a particular focus on endothelial cell damage by NETs. We comprehensively analyzed the differentially expressed genes (DEGs) in endothelial cells stimulated by NETs and took the intersection of these DEGs with MI-related genes. We expect to discover biomarkers of MI associated with NETs formation, especially regarding endothelial cell injury by NETs. Our findings may provide a new direction for studying NETs-mediated cell injury, particularly endothelial cells during MI.

Materials and Methods

Data Source

All gene expression data (microarray data, RNA-seq data) were downloaded from the Gene Expression Omnibus (GEO) database. The datasets GSE48060 and GSE66360 were annotated by platforms of GPL57, the dataset GSE62646 was annotated by platforms of GPL6244, GSE178883 and GSE179828 were annotated by platforms of GPL2467. In addition, we used an online web server (https://tabula-muris.ds.czbiohub.org/) to analyze gene expression in single-cell RNA-Seq data from a set of healthy mouse hearts.¹⁸

Identification and Enrichment Analyses of DEGs

GEO2R is a tool for differential analysis of expression profiling microarrays in the GEO database, with which we can compare two or more datasets of samples from the GEO database to obtain the DEGs. DEGs of GSE178883 and GSE179828 were respectively identified using GEO2R. The volcano diagram of DEGs was drawn using the “ggplot2” package. Based on the DAVID database (https://david.ncifcrf.gov/), the functional and molecular biological properties of the overlap were analyzed by GO (Gene Ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes) annotations.

Protein-Protein Interaction (PPI) Network Construction and Gene Set Enrichment Analysis (GSEA)

Based on the STRING database (http://string-db.org), the PPI network of overlap DEGs was created and visualized using Cytoscape software (version 3.9.1). GSEA was performed using the GSEA software (version 4.3.2). GSEA is a threshold-free method that examines all genes on the basis of their differential expression rank or other score, without first filtering the genes based on a predetermined set of criteria.

Evaluation of Immune Cell Infiltration

CIBERSORT, is a method for characterizing the cell composition of complex tissues from their gene expression profiles.¹⁹ The correlation analysis of 22 infiltrating immune cell types was evaluated using the “corrplot” package in the R software. ssGSEA is an extension of the GSEA method that allows for the definition of enrichment scores that represent the absolute enrichment of genomes in each sample in a given dataset.²⁰ xCell, based on the ssGSEA method, performed cell type analysis using gene expression data from 64 immune and stromal cells.²¹

Weighted Gene Co-Expression Network Analysis (WGCNA)

WGCNA is a computational approach in systems biology that designed to characterize the patterns of correlation between genes across several microarray datasets.²² WGCNA can be used to identify modules of highly correlated genes, to summarize these modules using the module eigengene or an intramodular hub gene, and to calculate module membership measures. In this study, the gene expression matrices were subjected to hierarchical cluster analysis in order to identify and remove any outliers. Subsequently, a correlation heatmap was generated to visualize the relationship between modules and traits, with the inclusion of P-values within the respective boxes. The relationship between gene expression and AMI was determined by evaluating gene significance and module membership.

Machine Learning

To find potential hub genes, three different machine learning algorithms were utilized. The Lasso regression analysis technique, performed through the “glmnet” R package, utilizes regularization to minimize prediction error. Support vector machine (SVM) is a supervised classification algorithm based on statistical learning theory. The random forest (RF) algorithm, implemented using the “randomForest” R package, was employed to identify a subset of candidate genes. This algorithm leverages ensemble learning by integrating numerous trees to enhance predictive accuracy. The overlapping genes identified by the Lasso model, SVM-RFE and random forest were defined as hub genes for subsequent research and validation.

Animals Experimental Model

All experimental animal procedures were approved by the Animal Care and Use Committee of Renmin Hospital of Wuhan University, and were also in accordance with the Guidelines for the Care and Use of Laboratory Animals published by the US National Institutes of Health. The mice were placed on the heating operation pad in the right lateral decubitus position, near the left axilla, the mouse skin and the third fourth rib were cut, the heart and left anterior descending artery (LAD) were fully exposed, the blood vessels were ligated with 7–0 stitched sutures at approximately 2 mm from the lower margin of the left atrial appendage under a stereomicroscope, and the distal area of the ligature was whitish due to ischemia, which was indicated by successful surgery. Mice in the sham group took the same operation but without suture ligation to ligate the blood vessels, and then the muscle and skin at the incision were sutured by layer.

Echocardiography and Speckle-Tracking Analysis

The day before mice were sacrificed, they were anesthetized with isoflurane (1.5–2%) and underwent transthoracic echocardiography to assess cardiac function. Echocardiographic parameters, such as left ventricular internal dimension diastolic (LVIDd), left ventricular internal dimension systolic (LVIDs), left ventricular ejection fraction (LVEF), heart rate (HR) and fractional shortening (FS), were measured using M-mode tracings. These measurements were obtained by averaging data from a minimum of five consecutive cardiac cycles. To evaluate the impact of myocardial remodeling following MI injury on cardiac dynamics at both global and regional levels, high-frequency echocardiography was performed using the speckle-tracking algorithm. To assess radial and longitudinal strain, strain rate, displacement, and velocity, cardiac cycles were acquired from the parasternal long-axis view.

Immunohistochemistry (IHC)

Isolated mouse hearts were immediately placed in 10% KCl solution and fixed in 10% formalin for 12 hours, followed by gradient alcohol dehydration, xylene hyalinization and paraffin embedding, and made into sections (4–5 μm). The paraffin sections were fully immersed in xylene solution, and after the immersion was completed, they were removed and dehydrated in a gradient alcohol solution. Paraffin sections were heated for antigen recovery using the pressure cooker method, and the sections were blocked for 1 h and then incubated overnight in a refrigerator at 4°C with antibodies against IL-1β (Abmart, TA5103), MMP9 (ABclonal, A11521), and the NFIL3 (Proteintech, 11773-1-AP) receptor, followed by incubation of the anti-rabbit/mouse EnVisionTM +/horseradish peroxidase reagent at a temperature of 37 °C for an additional 1 h, and then detected using a DAB staining kit. Images were quantified with ImageJ software.

RNA Isolation and Quantitative Real-Time PCR

Total RNA was extracted from the left ventricular tissue samples by grinding and lysis with TRIzol lysis solution, extracted by trichloromethane, precipitated by isopropanol and rinsed by 75% ethanol. The purity as well as the concentration of the proposed RNA was determined by NanoDrop One/Onec (Thermo Scientific), and the dilution volume was calculated and reverse transcription was performed. Light Cycler 480 SYBR Green 1 Master Mix was used to perform qRT-PCR. The data was analyzed using the 2^-(∆∆Ct) method and normalized to GAPDH. The details about all primer sequences are listed in Table S1.

Statistical Analysis

Statistical analyses were performed with GraphPad Prism software (version 9.0) or R software. The experimental data is in the form of mean ± SD. Two independent sample t-tests are used for sample comparison between two groups, and one-way ANOVA is used for sample comparison between multiple groups. A P value less than 0.05 was considered statistically significant.

Results

Immune Cell Infiltration Analysis

The research flowchart is presented in Figure 1. Cardiac healing following MI involves the mobilization and activation of immune cells.²³ For example, monocytes circulate the vasculature at a steady state and are recruited to sites of inflammation and engage in the development of cardiovascular disease.²⁴ Therefore, we performed an immune infiltration analysis of the gene expression dataset GSE48060 using two algorithms, “CIBERSORT” and “ssGSEA” to determine whether activated neutrophils are linked to the different disease states. The correlations between the infiltrating immune cells of “CIBERSORT” are shown in Figure 2A. The immune cell composition of samples from the AMI and control groups was compared (Figure 2B). Boxplot was generated to visualize the differences in immune infiltrating cell types of “CIBERSORT” between the AMI and control groups (Figure 2C). The boxplot analysis of the differential in immune cell infiltration demonstrated that patients with AMI exhibited a significantly elevated amounts of neutrophils compared to the control group. We also used the ssGSEA algorithm to analyze the immune cell infiltration in GSE48060 (Figure 2D and E). This algorithm compares the expression of 28 immune cell subtypes between the AMI and control groups. Similar to the results of the CIBERSORT analysis, there were significant differences in neutrophil expression levels between the two groups.

Figure 1 Flowchart of the study design. Abbreviations: AMI, acute myocardial infarction; PMA, phorbol myristate acetate; NETs, neutrophil extracellular traps; NRGs, NETs-related genes; DEGs, differentially expressed genes; HUVECs, human umbilical vein endothelial cells.

Figure 2 Immune cell infiltration analysis. (A) Correlation matrix of 20 immune cell subtype compositions. Blue represents positive correlation, Orange represents negative correlation, and the number in the square represents correlation. (B) Proportion of 20 subtypes of immune cells in different samples from AMI and control groups. (C) The proportion of the 20 immune cell subtypes in the AMI group and the control group is visualized in a barplot. (D) Heatmap of immune cell subtypes expression in AMI and control groups. (E) The barplot compares 28 immune cell subtypes proportion between AMI and control groups based on ssGSEA algorithm. *P < 0.05; **P < 0.01; ***P < 0.001.

Functional Annotation of PMA-Stimulated Neutrophils

NETs are reticular structures released by neutrophil activation and are closely related to the subsequent damage after MI.¹⁷ NETs are involved in the inflammatory pathophysiology of a variety of illnesses, such as cardiovascular and infectious diseases.²⁵ PMA can induce neutrophil aggregation and subsequent NETs release. Therefore, we analyzed differential gene expression after PMA stimulation of neutrophils. We selected human circulating neutrophils stimulated with 100 nM PMA or treatment with a control medium in GSE178883. Volcano plot showing the expression of differential genes between the two groups (Figure 3A). Figure 3B shows the top GO items under biological process (BP), molecular function (MF), and cellular component (CC). Next, functional enrichment analysis was performed based on the DEGs. The KEGG pathway enrichment analysis showed that the up-DEGs were mainly enriched in antigen processing and presentation, Th17 cell differentiation and Th1 and Th2 cell differentiation (Figure 3C). The down-DEGs were mainly enriched in neutrophil extracellular trap formation, Fc gamma R−mediated phagocytosis and cytokine−cytokine receptor interaction (Figure 3C). The results of GSEA analysis showed that HALLMARK_MYC_TARGETS_V2, HALLMARK_TNFA_SIGNALING_VIA_NFKB and HALLMARK_MYC_TARGETS_V1 were the most enriched gene sets in PMA group (Figure 3D). This confirmed that activation of inflammation-related pathways is an important feature of PMA-stimulated neutrophils.

Figure 3 Functional annotation of PMA-stimulated neutrophils. (A) Volcano plot of the DEGs in GSE178883. (B) GO enrichment analysis of DEGs. (C) KEGG functional enrichment analysis of up- and down-regulated genes in GSE178883. (D) GSEA comparing PMA treatment with vehicle control. Enrichment plots of key pathways in PMA treatment neutrophils. Abbreviations: PMA, phorbol myristate acetate. NES, Normalized Enrichment Score.

Significant Module Genes Identification in AMI via WGCNA

The WGCNA method was used to identify gene-gene modules associated with AMI and controls. We first eliminate the abnormal samples by setting a threshold (Figure 4A). A scale-free co-expression network was then established with the soft thresholding power to 14 (Figure 4B). A total of five modules were identified for future analysis. The primed modules were eventually displayed under the clustering tree (Figure 4C). We made a map of the connections between the identified modules, and the brown module showed a strong correlation with AMI (Figure 4D). The correlation between the module mean and gene significance (GS) for the brown module was found to be highly significant, providing further evidence suggesting that genes inside the brown module may possess functional significance related to AMI (Figure 4E). Figure 4F shows the heatmap of the selected gene. A total of 311 genes within the brown module were included in the subsequent study (Table S2).

Figure 4 Significant module genes identification in AMI via WGCNA. (A) Sample clustering to detect outliers (B) Mean connectivity analysis and scale-free index analysis for different soft-threshold power. (C) A dendrogram illustrating the clustering of genes, with distinct modules represented by different colors. (D) Heatmap of the correlation between module characteristic genes and AMI (E) A scatter plot was generated to show the association between GS and module membership (MM) within the brown module. (F) Network heatmap plot of genes.

Identification of Hub Genes Associated with AMI and NRGs Through Machine Learning

We next obtained NETs-related genes (NRGs) from the previous literature.²⁶ After taking the intersection of NRGs and MI-related genes from the brown module, we finally obtained 28 Intersecting genes (AMI-NRGs) (Figure 5A). The KEGG pathway enrichment analysis showed that IL-17 signaling pathway, TNF signaling pathway, Toll-like receptor signaling pathway and NF-kappa B signaling pathway were significantly enriched (Figure 5B). The results of GO functional enrichment indicated that the AMI-NRGs were mainly enriched in “inflammatory response”, “positive regulation of inflammatory response”, “neutrophil chemotaxis” regarding BP; “extracellular region”, “extracellular space”, “specific granule membrane” regarding CC; “identical protein binding”, “chemokine activity”, “RAGE receptor binding” regarding MF (Figure 5C). We then imported the AMI-NRGs into the STRING database and built a PPI network via Cytoscape (Figure 5D). Subsequently, three machine learning algorithms, including SVM-RFE, random forest, and Lasso regression, were employed to identify and select potential hub genes (Figure 5E–H), and we finally got four hub genes (FCAR, IL1B, MMP9 and NFIL3) (Figure 5I).

Figure 5 Identification of hub genes associated with AMI and NRGs through machine learning. (A) The 28 overlay genes between WGCNA analysis and NETs-related genes were presented by Venn diagram. (B) The most enriched KEGG pathways were presented by the Sankey diagram. (C) GO enrichment analysis of the AMI-NRGs. (D) PPI network of AMI-NRGs, isolated nodes were removed. (E and F) Visualization of Lasso regression. (G) The column shows that the random forest algorithm was used to order genes according to their importance score. (H) The hub genes were screened by the SVM-RFE algorithm. (I) The Venn diagram displays four common genes identified by Lasso, random forest and SVM-RFE algorithms, which were identified as the hub genes in AMI-NRGs.

Screening of DEGs from GSE179828

Neutrophil-mediated endothelial cell activation and injury play important roles in multiple disease states.²⁷ The release of NETs in blood vessels may trap red blood cells (RBCs), activate platelets, and damage the endothelium thereby promoting coagulation, vascular occlusion, and thrombosis.²⁷ Therefore, we analyzed the RNA-sequencing data of NETs-stimulated endothelial cells. We finally identified a set of 985 genes that were differentially expressed in human umbilical vein endothelial cells (HUVECs), 623 genes are upregulated, and 362 genes are downregulated (Figure 6A). The KEGG pathway enrichment analysis showed that TNF signaling pathway, NF-kappa B signaling pathway and Osteoclast differentiation were significantly enriched (Figure 6B). The GO analysis showed that the differential genes were most significantly enriched in the regulation of transcription from RNA polymerase II promoter in BP, nucleus in CC, and protein binding in MF (Figure 6C–E, Table S3).

Figure 6 Screening of DEGs from GSE179828. (A) Volcano plot of GSE179828. (B) The KEGG enrichment barplot of the DEGs. (C–E) The top 20 GO items under BP, CC, and MF.

Identification of Hub Genes Associated with AMI and NETs Stimulated Endothelial Cells

We next take the intersection of the endothelial cells-related genes acquired from GSE179828 with AMI-related genes from the brown module and finally obtained 59 intersecting genes (AMI-ECRGs) (Figure 7A). Functional enrichment analysis was performed based on the AMI-ECRGs (Figure 7B and C). The enrichment results are similar to the analysis outcome of AMI-NRGs, TNF signaling pathway, NF-kappa B signaling pathway and IL-17 signaling pathway are the top three enriched pathways (Figure 7B). We then imported the AMI-ECRGs into the STRING database and built a PPI network via Cytoscape (Figure 7D). The CytoHubba plugin’s four algorithms (Degree, Stress, Betweenness, and Radiality) were then utilized to identify common genes (Figure 7E). Finally, we screen out potential hub genes using three machine learning algorithms (random forest, Lasso regression, and SVM-RFE) (Figure 7F–I). CXCL2, ICAM1, and ZFP36 were identified as the intersection of common genes from SVM-RFE, random forest, and Lasso. (Figure 7J).

Figure 7 Identification of hub genes associated with AMI and NETs stimulated endothelial cells. (A) The 59 overlay genes between WGCNA analysis and the differentially expressed gene of GSE179828 were presented by Venn diagram. (B) The most enriched KEGG pathways were presented by the Sankey diagram. (C) GO enrichment analysis of the AMI-ECRGs. (D) PPI network of AMI-ECRGs, isolated nodes were removed. (E) The intersection of the first 15 genes from each of the four algorithms (Degree, Stress, Radiality and Betweenness) in the cytoHubba plugin resulted in the identification of eleven hub genes. (F and G) Visualization of Lasso regression. (H) The column shows that the random forest algorithm was used to order genes according to their importance score. (I) The hub genes were screened by the SVM-RFE algorithm. (J) The Venn diagram displaying three common genes identified by Lasso, random forest and SVM-RFE algorithms, which were identified as the hub genes in AMI-ECRGs.

Expression Profiles of Hub Genes in Different Cell Types of the Heart

The network of hub genes and their co-expression genes was created using the GeneMANIA platform, as shown in Figure 8A and B. Next, we used a publicly available single-cell RNA-Seq and FACS database of healthy mouse hearts to verify hub gene expression (Figure 8C and D). Because FCAR is not expressed in mice, we selected the remaining six hub genes for analysis. As shown in Figure 8D, Il1b, Cxcl2 and Mmp9 were primarily expressed in leukocytes and icam1 was mainly expressed in endothelial cells. Nfil3 is mainly expressed in fibroblasts and cardiomyocytes, whereas Zfp36 is ubiquitously expressed in a wide range of cells including endothelial cells, fibroblasts, and leukocytes. The dataset GSE95755 serves as a comprehensive transcription resource that encompasses diverse cardiac cell types throughout the processes of heart development, repair, and regeneration.²⁸ We selected cells isolated from infarcted and non-infarcted areas of adult mouse hearts for our subsequent validation. Figure 8E shows the expression levels of hub genes in the four cell types (cardiomyocyte, fibroblast, endothelial cell and leukocyte) between the MI and sham groups. The expression levels of Il1b, Mmp9, Cxcl2 and Zfp36 were elevated in all four cell types after the MI operation compared to the sham group (Figure 8F–I).

Figure 8 Expression profiles of hub genes in different cell types of the heart. (A) Hub genes include AMI-NRGs (FCAR, IL1B, MMP9 and NFIL3) and AMI-ECRGs (CXCL2, ICAM1 and ZFP36). (B) Hub genes and their co-expression genes were analyzed by GeneMANIA. (C and D) Expression of hub genes in different cell types of the healthy mouse heart. Measured by FACS in the Tabula Muris publically-available dataset. (E) Expression levels of hub genes in different cell types between MI and sham groups. (F–I) Expression levels of hub genes in (F) cardiomyocyte, (G) fibroblast, (H) leukocyte and (I) endothelial cells between MI and sham groups.

The Diagnostic Value of Hub Genes

We evaluated the predictive accuracy of the hub genes (AMI-NRGs and AMI-ECRGs) in distinguishing normal individuals an AMI. Receiver operating curve (ROC) were used to assess the diagnostic value of each candidate hub gene. The normal individuals and MI patients in the dataset of GSE66360 could be accurately distinguished (Figure 9A–D and K–N). The AMI-NRGs and AMI-ECRGs also demonstrated significant accuracy in distinguishing AMI from normal individuals in GSE66360, with an AUC of 0.938 (95% CI 0.8911 to 0.9848) and an AUC of 0.909 (95% CI 0.8547to 0.9633) (Figure 9E and N). We also evaluated our hub genes in another dataset GSE62646 (Figure 9F–J and O–R) and the AMI-ECRGs got good results (Figure 9R). Additionally, we use several datasets including GSE66360, GSE62646 and GSE48060 to validate hub gene expression (Figure 9S). The relative expression levels of CXCL2, FCAR, ICAM1, MMP9, and NFIL3 were significantly higher in AMI than in the control groups in GSE66360 and GSE48060.

Figure 9 The diagnostic value of hub genes. (A–D and F–I) The ROC curve of FCAR, IL1B, MMP9 and NFIL3 in (A–D) GSE66360 and (F–I) GSE62646. (E and J) The ROC curve of AMI-NRG in (E) GSE66360 and (J) GSE62646. (K–M and O–Q) The ROC curve of CXCL2, ICAM1 and ZFP36 in (K–M) GSE66360 and (O–Q) GSE62646. (N and R) The ROC curve of AMI-ECRG in (N) GSE66360 and (R) GSE62646. (S) The expression of hub genes in GSE66360, GSE62646 and GSE48060. The comparison between the AMI and control groups used the mean t-test. P < 0.05 was considered to be statistically significant. *P < 0.05, **P < 0.01, ***P < 0.001.

Correlation Analysis Between Hub Genes and Infiltrating Immune Cell Types

To conduct a more comprehensive evaluation of the differences in immune cell infiltration and signature gene sets between AMI and control samples, the CIBERSORT, ssGSEA and xCell algorithms were employed. Based on the correlation study conducted using CIBERSORT, it was found that, except CXCL2, the remaining six hub genes exhibited a significant positive correlation with the infiltration of neutrophils and macrophages M0 (Figure 10A). Consistent with the results of the analysis using CIBERSORT, the ssGSEA results showed that except CXCL2 the other six hub genes were positively correlated with activated dendritic cell, eosinophil, macrophage, mast cell, neutrophil and plasmacytoid dendritic cell (Figure 10B). In the test dataset GSE62646, our analysis also showed that neutrophils and monocytes were positively correlated with hub genes (Figure 10C and D). Next, we analyzed the dataset using the xCell algorithm. The results showed a high positive correlation of genes with monocytes, NKT and neutrophils (Figure 10E and F).

Figure 10 Correlation analysis between hub genes and infiltrating immune cell types. (A, B) Correlation heatmap showing the relationship between seven hub genes and the infiltrated immune cells based on CIBERSORT and ssGSEA algorithms in GSE48060. (C, D) Correlation heatmap showing the relationship between seven hub genes and the infiltrated immune cells based on CIBERSORT and ssGSEA algorithms in GSE62646. (E, F) Correlation heatmap showing the association between seven hub genes and the infiltrated immune cells based on the xCell algorithm. *P < 0.05, **P < 0.01, ***P < 0.001.

Establishing a Model of MI

To verify the expression level of hub genes, we established a MI model using C57BL/6J mice. Echocardiography was conducted to evaluate the functional alterations of the mice in each experimental group. Figure 11A shows the representative B- and M-mode echocardiographic imaging of the heart. Mice after MI surgery exhibited significantly decreased LVEF and FS but exhibited significantly increased LVIDs and LVIDd compared with the sham controls (Figure 11B–F). Strain analysis and 2-dimensional speckle-tracking echocardiography were used to further evaluate the detrimental cardiac remodeling that occurs in mice following an MI. As shown in Figure 11G–J, an evaluation of the function of the left ventricular using longitudinal and radial strain analysis showed a notable decrease in global longitudinal strain (GLS) and lower radial strain in the anterior segments in post-MI hearts of mice compared to mice in the sham group.

Figure 11 Establishing a model of myocardial infarction. (A) Representative B- and M-mode echocardiographic images of sham and MI mice. (B–F) Echocardiographic quantification of LVEF, FS, HR, LVIDd and LVIDs in the indicated groups (n=6). (G) Representative speckle-tracking echocardiography-based analysis of left ventricular function in sham and MI mice at day 7 post-MI. (H) Average radial strain was measured by spot tracking technology (n=5). (I) Representative long-axis radial myocardial strain maps deconvoluting infarction-induced regional myocardial wall deformation in sham and MI mice at day 7 post-MI. (J) GLS was measured by spot tracking technology (n=5). ***P < 0.001, ****P < 0.0001.

Validation of Gene Expression

To further demonstrate the expression pattern of hub genes, qRT-PCR verification was performed in the left ventricular tissue samples of the sham and MI groups. The results indicated that IL-1β, MMP9 and NFIL3 were highly expressed in MI samples (Figure 12A and B). However, the expression levels of ICAM, ZFP36, and CXCL2 were not significantly different between the two groups, which may be due to our insufficient sample size and single time point validation after MI. Next, we used immunohistochemistry to further verify the gene expression, in parallel, IL-1β, MMP9 and NFIL3 showed substantial upregulation in the infarcted areas (Figure 12C). Statistical analysis showed a significant increase in positive areas in the infarcted area, but no significant difference between the non-infarcted area and the sham groups (Figure 12D-F).

Figure 12 Validating gene expression. (A) The mRNA expression of AMI-NRGs (IL-1β, MMP9, NFIL3) was identified via qRT-PCR (n = 5–6). (B) The mRNA expression of AMI-ECRGs (ICAM1, ZFP36, CXCL2) were identified via qRT-PCR (n = 5–6). (C) Representative images of IHC staining of hearts from MI and sham groups mice. (D-F) Statistical analysis of IHC staining (n =3). Scale bar, 100μm. *P < 0.05, **P < 0.01, ***P < 0.001.

Discussion

There is increasing evidence that supports the significant involvement of neutrophils in the progression of coronary artery disease (CAD), particularly in the subsequent problems that arise, such as acute coronary syndrome and heart failure.²⁹ For example, activated lesional smooth muscle cells (SMCs) attract neutrophils, triggering the ejection of neutrophil NETs.³⁰ Among them, the nuclear protein histone H4 binds to and lyses SMCs, this process can ultimately lead to the destabilization of plaques.³⁰ Therefore, neutrophils may be an independent predictor for risk stratification in patients with ACS and/or cardiac revascularization when analyzed concomitantly with other markers of inflammation.³¹

In the present study, we first analyzed the immune infiltration of MI based on the CIBERSORT and ssGSEA algorithms. Similar to the results of previous studies,³² we found that neutrophils were higher in the AMI group. In a subsequent enrichment pathway analysis of the DEGs in GSE178883, we found that in addition to the formation of NETs, the down-regulated DEGs were also enriched in Fc gamma R−mediated phagocytosis and Cytokine−cytokine receptor interaction. This may reflect neutrophils as one of the major players during acute inflammation.³³ We also observed that up-regulated DEGs were enriched in Th17 cell differentiation and Th1 and Th2 cell differentiation. This finding suggests that, apart from their inherent function in resisting infections and responding to injuries, neutrophils are increasingly recognized as crucial regulators of the adaptive immune system.³⁴

Using three machine learning algorithms, we finally filtered out four hub genes in the AMI-NRGs (FCAR, IL1B, MMP9 and NFIL3). The FCAR gene has been identified as a member of the immunoglobulin gene superfamily and is responsible for encoding a receptor that specifically binds to the Fc region of immunoglobulin A (IgA). FcαRI/CD89 is unique in that it is neither expressed in mice nor is there any homolog in mice.³⁵ The previous studies have shown that neutrophils possess natural anti-tumor properties and can induce potent tumor cell killing via targeting the immunoglobulin A Fc receptor.³⁶ Here, our study suggests that neutrophils may mediate cellular injury during MI through the immunoglobulin A Fc receptor. Similar to FCAR, IL-1β, recognized as a member of a much larger IL-1 family of cytokines, is expressed primarily in myeloid cells.³⁷ Large numbers of neutrophils can quickly accumulate at infection or tissue damage sites, and become the predominant local source of IL-1.³⁸ IL-1β can be converted from an inactive precursor to a biologically active cytokine through inflammasome activation of the intracellular cysteine protease caspase-1 and caspase-1 independent cleavage mediated by serine proteases derived from neutrophils.³⁹ Many studies have demonstrated the important role of IL-1 in MI.⁴⁰ Thus, inhibition of IL-1β may reduce inflammatory injury in the acute phase of MI by inducing an anti-inflammatory effect.

In neutrophils, preformed MMP-9 is mainly stored in gelatinase granules and subsequently released in response to stimulation by inflammation or tissue injury.⁴¹ In the early stages after MI, MMP9 is mainly derived from neutrophils.⁴² Serum exosome MMP9 expression is elevated in patients with MI, and can be an effective biomarker for diagnosis of AMI.⁴³ Nishiguchi et al demonstrated that local MMP-9 level could determine the early clinical presentation in patients with acute MI.⁴⁴ Our experiments also confirmed that MMP9 is highly expressed in MI tissues. NFIL3, also known as Adenovirus E4 promoter-binding protein (E4BP4), is one of the mammalian basic leucine zipper transcription factors. Research has revealed that NFIL3 is linked to numerous physiological and biochemical processes.⁴⁵ For example, NFIL3 has been found to be related to immune-mediated diseases, like SLE, arthritis and Crohn’s disease.⁴⁶ In addition, some studies have demonstrated that NFIL3 may regulate the pathological process of heart failure through calcium signaling mechanisms, autocrine signaling, and insulin-like growth factor II receptor.⁴⁷ Our experiments confirm that NFIL3 is highly expressed in ischemic myocardial tissue and may be involved in neutrophil-mediated cell damage. Next, further studies are needed to determine the specific mechanism of NFIL3 in NETs-mediated cell injury.

To investigate possible targets of NETs-mediated endothelial cell injury during MI, we took the intersection of the brown module with the DEGs of GSE179828, and finally get three hub genes (CXCL2, ICAM1 and ZFP36). Although our experimental validation did not reveal significant differences in the expression levels of these genes in the control and MI groups, previous studies have shown the importance of these genes in the development of many cardiovascular diseases. For example, several studies found that CXCL2 has an important role in cardiovascular disease, especially in the development of atherosclerosis. The study conducted by Chang et al provided evidence supporting the notion that the absence of Sirt4 contributes to the progression of atherosclerosis through the activation of the NF-κB/IκB/CXCL2/3 signaling pathway.⁴⁸ Another study reveals that the joint suppression of TNFα, CXCL2, and CCL2 inhibits the ability of aged visceral fat transplants to promote the development of atherosclerosis.⁴⁹ The expression of ICAM-1 is increased in the milieu of atherosclerotic plaques, where it is observed to be expressed by endothelial cells, macrophages, and smooth muscle cells. Gross et al demonstrate that soluble intercellular adhesion molecule 1 (sICAM-1) concentration may be an early biomarker that indicates changes in the artery wall that accompany atherosclerosis.⁵⁰ ZFP36 is a small group of mRNA binding and destabilizing proteins expressed in almost all eukaryotes.⁵¹ A previous bioinformatics analysis study suggests that ZFP36 may be an endothelial cell senescence-related gene in patients with MI.⁵² Li et al also found significantly elevated ZFP36 protein levels in ischemic human heart tissue and knockdown of ZFP36 in HCMECs using siRNA significantly inhibited cell proliferation.⁵³ Our analysis suggests that ZFP36 may be one of the targets of endothelial cell injury by NETs during MI and may be a promising novel target for the treatment of ischemic heart disease.

Our study has several limitations. First, we did not integrate more datasets to screen for the DEGs, and more advanced bioinformatics and machine learning methods may screen for hub genes more accurately. Second, clinical samples and data are needed to further validate gene expression levels and their diagnostic value. Finally, further in vivo and in vitro loss-and-gain-of-function experiments are necessary to investigate the mechanisms of hub genes in NETs-mediated cellular damage.

Conclusions

In conclusion, our study identified seven candidate diagnostic hub genes, including three genes associated with NETs-related endothelial cell injury. These hub genes are likely to serve as important markers for the diagnosis of MI. Furthermore, analyzing these hub genes and immune cells may provide new insights into potential therapeutic approaches for NETs-related cell injury during MI.