Identification of diagnostic immune-related gene biomarkers for predicting heart failure after acute myocardial infarction

Abstract Post-myocardial infarction heart failure (HF) is a major public health concern. Previous studies have reported the critical role of immune response in HF pathogenesis. However, limited studies have reported predictive immune-associated biomarkers for HF. So we attempted to identify potential immune-related indicators for HF early diagnosis and therapy guidance. This study identified two potential immune-related hub genes (IRHGs), namely CXCR5 and FOS, using bioinformatic approaches. The expression levels of CXCR5 and FOS and their ability to predict long-term HF were analyzed. Functional enrichment analysis revealed that the hub genes were enriched in immune system processes, including the interleukin-17 and nuclear factor-kappa B signaling pathways, which are involved in the pathogenesis of HF. Quantitative real-time polymerase chain reaction revealed that the Fos mRNA levels, but not the Cxcr5 mRNA levels, were downregulated in the mice of the HF group. This study successfully identified two IRHGs that were significantly and differentially expressed in the HF group and could predict long-term HF, providing novel insights for future studies on HF and developing novel therapeutic targets for HF.


Introduction
Inflammation is involved in the pathogenesis of heart failure (HF).Under stress conditions, cardiomyocytes release proinflammatory cytokines [1], which trigger immune responses (including macrophage infiltration into the heart [2]) and mediate disease progression [3].T cells are reported to mediate the progression of pressure overload-induced HF [4][5][6].Additionally, the immune response in the impaired myocardium is characterized by the activation of the classical human immune response, which is similar to the response to autoimmunity and infection [7].These findings indicate the potential correlation between myocardial diseases, cytokines, and immune cell populations.Previous studies have demonstrated the role of immune response in HF pathogenesis [8].However, limited studies have identified the immune-associated biomarkers for predicting HF.So we attempted to identify potential immune-related indicators for HF early diagnosis and therapy guidance.
Recently, bioinformatics approaches, a type of machine learning, have been used to identify biomarkers for the diagnosis and therapy guidance of HF [9,10].Cell-type identification by estimating relative subsets of RNA transcripts (CIBERSORT), a widely used algorithm for identifying immune cells, enables the assessment of relevant subpopulations of RNA transcripts and the quantification of cell components [11].Additionally, CIBERSORT enables the elucidation of the immune cell landscape in tumors [12,13].
Bioinformatics approaches have identified several candidate genes for post-acute myocardial infarction (MI) (AMI) HF.However, the immune mechanisms of these genes in post-AMI HF have not been elucidated.Maciejak et al. [14] revealed potential biomarkers for predicting post-AMI HF, including FMN1, JDP2, and RNASE1.Similarly, Li et al. [15] discovered candidate time series differentially expressed genes (DEGs) associated with post-AMI HF, namely FADS2, LRRN3, GPR15, and AK5.
This study aimed to identify immune-associated biomarkers to complement the known predictive indicators [14,15] and reveal the differential immune cell infiltration statuses between the HF and non-HF groups.Additionally, this study explored the diagnostic value of the identified biomarkers to provide novel insights for future studies on HF and develop novel therapeutic targets for HF.

Data acquisition
The microarray expression dataset GSE59867 was downloaded from the Gene Expression Omnibus (www.ncbi.nlm.nih.gov/geo) database using the "BiocManager" package in R. The GSE59867 dataset was developed based on the GPL6244 platform ([HuGene-1_0-st] Affymetrix Human Gene 1.0 ST Array [transcript (gene) version]).The dataset comprised the data of 111 patients with ST-elevation MI at four different time points (admission, the day of AMI diagnosis; discharge, days 4-6 post-AMI; month 1 post-AMI; and month 6 post-AMI).The control group comprised 46 patients who were diagnosed with stable coronary artery disease on the day of admission.After 6 months, 9 patients were considered to exhibit HF and 8 patients did not exhibit HF according to the first and fourth quartiles of plasma N-terminal pro-B type natriuretic peptide (NT-proBNP) level and left ventricular ejection fraction.The demographic characteristics, including age, sex, body mass index, medical history, drug administration, and smoking history, were not significantly different between the HF and non-HF groups (p > 0.05).Immune-related genes (IRGs) were downloaded from the ImmPort database (https:// www.immport.org/home).

Data pre-processing
The GSE59867 raw dataset was subjected to pre-processing methods, including background calibration, normalization, and log 2 transformation, using the "Affy" package [16] in R "Bioconductor."Based on the probe annotation information, the probes were converted into gene symbols.If multiple probes corresponded to the same gene, the median value was considered its gene expression value.Additionally, the probes matching with more than one gene were removed.

DEG analysis
DEGs between the HF and non-HF groups were screened using the "Limma" package [17] in R based on the following criteria: p < 0.05 and |log 2 -fold change (FC)| > 0.5.DEGs with log 2 FC > 0.5 were considered upregulated genes, whereas those with log 2 FC < −0.5 were considered downregulated genes.

Gene ontology (GO) and Kyoto
Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses and gene set enrichment analysis (GSEA) The potential functions of DEGs were evaluated using GO and KEGG pathway enrichment analyses with the "ClusterProfiler" package [18] in R language.The enrichment of GO terms and KEGG pathways was considered significant at p < 0.05.Functional enrichment of DEGs was further examined using GSEA as previously reported [19].GSEA was performed using the "ClusterProfiler" package in R language and the enrichment was considered significant based on the following criteria: |normalized enrichment score| > 1; p < 0.05; false discovery rate < 0.25.Gene symbols were converted to Entrez ID using the human genome annotation package "org.Hs.eg.db." and plotted using the "ggplot2" package [20] in R.

Protein-protein interaction (PPI) network construction and module analysis
STRING software (version 11.5) was used to construct the network of DEGs (|log 2 FC| > 0.5 and p < 0.05) with the combined score of >0.4 considered the cut-off criterion.
The PPI network was imported into Cytoscape software (version 3.8.2) for further visualization.The Molecular Complex Detection (MCODE) plugin (version 2.0.0) in Cytoscape software was used to identify the most significant module.The operating parameters were as follows: degree cut-off = 2; node score cut-off = 0.2; k-score = 2; and max.depth = 100.The top 10 hub genes were identified using the CytoHubba plugin (version 0.1) in Cytoscape.

Estimation of infiltrating immune cell abundance
The CIBERSORT algorithm enables the identification of various immune cells based on their expression profiles [11].
The normalized GSE59867 expression matrix was analyzed using the CIBERSORT algorithm to obtain the proportions of various immune cell types.The core algorithm of CIBE-RSORT is the deconvolution algorithm, which utilizes the single-cell gene expression profiles to compute the proportion of different cell types.

Receiver operating characteristic (ROC) curve analysis
The sensitivity and specificity of the hub genes to distinctly identify post-AMI HF were determined using ROC curve analysis.The pROC package in R language was used to obtain the sensitivity, specificity, and area under the ROC curve (AUC) values.

Quantitative real-time polymerase chain reaction (qRT-PCR)
Total mRNA was extracted using TRIzol reagent and reversetranscribed into complementary DNA using the Evo Maloney murine leukemia virus reverse transcription Premix for qPCR (Accurate Biotechnology, Hunan, China).qRT-PCR analysis was performed using SYBR ® Green Pro Taq HS Premix (Accurate Biotechnology, Hunan, China) and LightCysler480 Real-Time PCR System (Germany).The primers for Fos, Cxcr5, and Actb were synthesized by Tsingke Biology (Beijing, China).The gene expression levels were analyzed using the 2 −ΔΔCt method.The expression levels of target genes were normalized to those of Actb.The primers used in qRT-PCR analysis are listed in Table 1.

Western blotting
The samples were lysed using a mixture of radioimmunoprecipitation assay lysis buffer, phosphatase inhibitor, and phenylmethylsulfonyl fluoride (100:1:1) and an Ultrasonic Cell Disruptor.The lysates were subjected to electrophoresis.The resolved proteins were transferred to a polyvinyl difluoride membrane (Millipore, USA).The membrane was blocked with 5% milk powder in Tris-buffered saline containing Tween-20 (TBST) for 1 h.Next, the membrane was probed with the anti-ANP (1:1,000, Proteintech, China) antibodies at 4°C overnight, followed by incubation with the horseradish peroxidase-labeled secondary antibodies at room temperature for 1 h.After washing the membrane thrice with TBST, immunoreactive signals were developed using an enhanced chemiluminescence reagent (Fude Biological, Hangzhou, China).

Animal model establishment
C57BL/6 mice (22-25 g) obtained from the Animal Center of Nanfang Medical University were randomly divided into the sham (n = 5) and MI (n = 5) groups.Mice in the sham group underwent sham surgery, while those in the MI group were subjected to left anterior descending artery permanent ligation as described previously [21].Animals were maintained under the following conditions for 2 months: circadian cycle, 12 h light/dark cycle; temperature, 22°C; access to water and food, ad libitum.Mice were then anesthetized and sacrificed.

Masson's staining
The paraffin heart sections were deparaffinized, rehydrated, and subjected to Masson's staining with the Masson's kit (# Immune-related gene biomarkers in heart failure  3 G1340; Solarbio), following the manufacturer's instructions.The sections were sealed and observed.

Statistical analysis
All statistical analyses were performed using GraphPad Prism (version 8.0.2) and R software (version 4.2.4).Data are presented as mean ± standard error of the mean.The sensitivity and specificity of hub genes to distinguish post-AMI HF from post-AMI non-HF were calculated and represented using an ROC curve.All tests in this study were two-tailed.Differences were considered significant at p < 0.05.
Ethics approval and consent to participate: Animal experiments were performed according to the instructions of the Animal Ethics Committee of Nanfang Medical University and the ARRIVE guidelines.Efforts were made to minimize animal suffering.

Identification of DEGs between the non-HF and HF groups
To identify the critical genes involved in the pathogenesis of post-AMI HF, the GSE59867 dataset was analyzed, which revealed nine patients with HF and eight non-HF cases whose samples were collected at admission.The gene expression values in the GSE59867 dataset were normalized (Figure 1a and b).In total, 200 DEGs were identified between the non-HF and HF groups (84 upregulated genes  Immune-related gene biomarkers in heart failure  5 and 116 downregulated genes) (Figure 1c).The heatmap of DEGs was generated, which revealed a distinct expression pattern between the non-HF and HF groups (Figure 1d).To further assess the differential gene expression patterns between the HF and non-HF groups, the GSE59867 dataset was subjected to principal component analysis.The two groups were distinctly separated in the Dim1 component (Figure 1e).These findings indicate the well-replicated data between the HF and non-HF groups.

Functional enrichment analysis of DEGs
To examine the functions of the 200 identified DEGs, the DEGs were subjected to GO and KEGG pathway enrichment analyses using the "ClusterProfiler" packages in R language.The DEGs were mainly enriched in various GO terms as follows: biological process (BP) term (Figure 2a), immune system process; cellular component (CC) term (Figure 2b), external side of plasma membrane; molecular function (MF) term (Figure 2c), immune receptor activity, cytokine receptor activity, and cytokine binding.KEGG pathway (Figure 2d) enrichment analysis revealed that the DEGs were enriched in hematopoietic cell lineage.GSEA revealed that the top 10 enriched pathways in which the DEGs were enriched included amoebiasis and pantothenate and CoA biosynthesis (Figure 2e and f).

Potential key IRGs associated with HF
To examine the interaction between the identified DEGs, a PPI network was constructed using the STRING online tool (Figure 3a).The PPI network was imported into the Cytoscape software for further analysis.Module analysis was performed using the MCODE plugin to obtain the first module, which was subjected to subsequent analysis (Figure 3b).The first module was enriched in lymphocyte proliferation (Figure 3c).To identify the hub genes among the DEGs, the CytoHubba plugin in Cytoscape software was used to extract the top 10 hub genes (IL1B, CD28, CXCL8, IL2RA, KLRC4-KLRK1, CXCR5, CD40LG, FOS, TIMP1, and CCR6) (Figure 3d1 and d2).Functional enrichment analyses revealed that the DEGs were enriched in immune responses.Hence, the IRGs were downloaded from the ImmPort database.The intersection of 10 hub genes and IRGs revealed seven immune-related hub genes (IRHGs) (Figure 3e).

Functional enrichment analysis of IRHGs
IRHGs were subjected to GO and KEGG functional enrichment analyses.The IRHGs were enriched in different GO terms as follows: BP term, T-cell proliferation and lymphocyte proliferation (Figure 4a); CC term, external side of plasma membrane (Figure 4b); and MF term, cytokine activity (Figure 4c).KEGG pathway enrichment analysis revealed that IRHGs were mainly enriched in cytokine-cytokine receptor interaction, interleukin (IL)-17 signaling pathway, and nuclear factorkappa B (NF-κB) signaling pathway (Figure 4d).

Profiling of immune cells associated with HF
Next, the abundance of the immune cells in the AMI samples was profiled using the CIBERSORT algorithm (Figure 5a1 and a2).
Comparative analysis of the abundance of immune cells Immune-related gene biomarkers in heart failure  7 between the HF and non-HF groups revealed that the proportion of resting natural killer (NK) cells was significantly upregulated in the HF group (Figure 5b).

Validation of hub genes associated with HF
To further assess the potential values of the seven IRHGs, the expression levels of IRHGs and their ability to predict long-term HF were determined.CXCR5 and FOS were differentially expressed between the HF and non-HF groups (Figure 6a-g).Additionally, CXCR5 (AUC = 0.967; 95% confidence interval = 0.874-1) and FOS (AUC = 0.667; 95% confidence interval = 0.234-1) could predict longterm HF (Figure 6h and i).

Mouse experiments
The heart size in the MI group was higher than that in the sham group (Figure 7a).The representative images of heart sections from the sham and MI groups subjected to Masson's staining on day 60 post-surgery (scale bars = 500 μm) are presented in Figures 7b and c.Cardiac fibrosis was observed in the MI group (Figure 7c).The cardiac function of mice in the sham and MI groups was evaluated using Mmode echocardiography.The left ventricular movements in the MI group were significantly poor when compared with those in the sham group (Figures 7d and e).
The effects of HF post-AMI on the mRNA expression levels of Cxcr5, Fos, and atrial natriuretic factor (Nppa) were examined using qRT-PCR analysis.The Fos and Nppa mRNA levels in the MI group were upregulated when compared with those in the sham group.Meanwhile, the Cxcr5 mRNA levels were downregulated in the MI group (Figure 8a-c) (p < 0.05).The heart weight (HW), lung weight (LW), and tibia length (TL) were measured in each mouse.The ratios of HW/TL and LW/ TL were markedly different between the MI and sham groups (Figure 8d and e) (p < 0.05).Western blotting analysis revealed that the Nppa protein levels in the MI group were significantly higher than those in the sham group (Figure 8f and g) (p < 0.05).These results indicated that the genes identified in this study were associated with the development of HF post-AMI.This study revealed that CXCR5 and FOS were differentially expressed between the HF and non-HF groups and could predict long-term HF.Additionally, functional enrichment analysis revealed that CXCR5 and FOS were enriched in immune system processes associated with the progression of HF, including the IL17 and NF-κB signaling pathways.These novel biomarkers can complement current biomarkers for predicting long-term HF in patients with AMI and provide novel insights for understanding HF development.Previous studies have identified potential biomarkers for HF post-AMI, such as FMN1, JDP2, LRRN3, FADS2, GPR15, and AK5 [9,10].However, the correlation of these indicators with immune mechanisms has not been elucidated.This study revealed that the proportion of resting NK cells varied between the HF post-AMI and non-HF groups.Additionally, CXCR5 and FOS were demonstrated to predict long-term HF.Furthermore, in vivo experiments were performed to validate the in silico findings.Thus, this study demonstrated the role of immune response in the development of HF and provided a scientific rationale for future studies on HF.
CXCR5, a chemokine receptor of CXCL13, is expressed on mature B cells and follicular T helper cells and regulates homeostatic lymphoid cell trafficking and homing to B-cell follicles within secondary lymphoid organs [22,23].In combination with the members of the lymphotoxin/tumor necrosis factor family, CXCR5 and its ligand (CXCL13) promote the development and maintenance of secondary lymphoid organs.Cxcr5-deficient mice are reported to exhibit impaired B-cell follicle development [23,24].Heinrichs et al. [25] reported that compared with those in wild-type mice, the post-MI infiltration levels of B cells were downregulated Cxcr5-deficient mice or in mice treated with anti-CXCL13 neutralizing antibody.Recent studies have demonstrated that B cells and antibodies can modulate inflammation and remodeling after MI [25].Waehre et al. revealed that after 3 weeks of aortic banding, Cxcr5 knockout (Cxcr5 −/− ) mice exhibited left ventricular dilation when compared with wild-type mice.Additionally, the myocardial CXCR5 levels are upregulated in patients with HF [26].One study examined CXCR5 expression in thrombus obtained from plaque rupture during MI and reported strong CXCR5 immunostaining of CXCR5 in the thrombus of patients with ST-elevation MI and unstable carotid lesion.Additionally, factors released from thrombin-activated platelets upregulated the levels of CXCR5 in monocytes [27].At day 7 post-MI, MMP12 inhibitor upregulated CXCR5 [28].Thus, the activation of immune responses is critical for the pathogenesis of MI and HF, such as promoting adverse cardiac remodeling and left ventricular disorder.
Palomer et al. reported that miR-146a overexpression downregulated the FOS mRNA and protein levels.MMP9 activity was suppressed upon miR-146a-mediated downregulation of the FOS-AP-1 pathway.HF is closely associated with the upregulation of MMP9 protein [29].Xue et al. reported that the anti-cardiac fibrotic effect of miR-29b-3p is partly mediated by targeting FOS [30].Thus, FOS is involved in the pathogenesis of cardiac fibrosis and HF.
Immune cells are reported to promote myocardial injury [31][32][33].In particular, monocytes and macrophages in the infarcted heart can elicit inflammatory responses and promote left ventricular adverse remodeling and HF [32,34,35].This study demonstrated that the proportions of resting NK cells were significantly upregulated in patients with HF.These findings improved our understanding of the novel mechanism underlying post-AMI HF.
KEGG pathway enrichment analysis revealed that IRHGs were enriched in the IL17 and NF-κB signaling pathways.The IL17 pathway, which is a potential therapeutic target for immune-mediated diseases, has been widely studied and demonstrated to be involved in the pathogenesis of cardiac diseases [36,37].Compared with those in control subjects, the IL17 levels are upregulated in patients with HF [36,38].IL17 promotes myocardium apoptosis and ischemic myocardium remodeling [39].A previous study reported that IL17 promotes fibrosis, collagen generation, and apoptosis response and that these effects of IL17 were reversed upon treatment with neutralizing antibodies [37].Il17 is upregulated during the pathogenesis of HF in mice through the NF-κB-dependent suppression of Atp2a2 and Cacna1c expression [36].
The transcription factor NF-κB is the downstream mediator of the IL17 signaling pathway [40].The induction of NF-κB promotes the release of inflammatory cytokines, chemokines, and adhesion molecules from innate and adaptive immune cells.The subsequent inflammatory cascades adversely affect the cardiovascular system [41].NF-κB, which mediates the pathological processes of various cardiovascular diseases, is reported to be involved in promoting inflammation, cell survival, cell differentiation, and cell proliferation [41], which contribute to the progression of cardiovascular diseases, such as HF, myocardial ischemia, hypertension, and atherosclerosis [41].
This study has several limitations.The sample size in this study was small, which may lead to biased conclusions owing to population differences.The CIBERSORT algorithm has several limitations.For example, the CIBERSORT algorithm can systematically underestimate or overestimate certain cell types, even though it is associated with a relatively lower estimation bias than other similar methods [16].
Our study successfully identified two IRHGs that were differentially expressed between the HF and non-HF groups.These genes could predict long-term HF.Thus, the findings of this study provided novel insights for understanding HF development.

Figure 1 :
Figure 1: (a) The expression values before normalization.(b) The expression values after normalization.(c) The volcano map of DEGs (red and blue dots represent upregulated and downregulated genes, respectively).(d) The heatmap of the DEGs (each row represents one DEG, while each column represents one sample.The red and blue colors indicate upregulated and downregulated DEGs, respectively).(e) Principal component analysis of DEGs between the HF and non-HF groups revealed that the two groups were separated in the Dim1 component.

Figure 2 :
Figure 2: GO and KEGG pathway enrichment analyses and GSEA of DEGs.The red and blue colors represent high and low adjusted p values, respectively.The size of the dot indicates the amount of enriched DEGs.(a) The top 10 BP enrichment terms.(b) The top 10 CC enrichment terms.(c) The top 10 MF enrichment terms.(d) Enriched KEGG pathways.(e and f) The top 10 enriched pathways were identified using GSEA.

Figure 3 :
Figure 3: (a) The PPI network of DEGs was constructed using the STRING online tool.(b) The first module obtained using the Cytoscape software and MCODE plugin.(c) Enrichment analysis of the first module.(d1 and d2) The top 10 hub genes among DEGs were identified using the CytoHubba plugin in Cytoscape software.(e) Venn diagram shows overlapping genes between hub genes and IRGs.

Figure 4 :
Figure 4: GO and KEGG analyses of IRHGs.The red and blue colors represent high and low adjusted p values, respectively.The size of the dot indicates the amount of enriched IRHGs.(a) The top 10 BP enrichment terms.(b) The top 10 CC enrichment terms.(c) The top 10 MF enrichment terms.(d) The enriched KEGG pathways.

Figure 5 :
Figure 5: (a) The proportions of immune cells in the HF (a1) and non-HF groups (a2).(b) Comparative analysis of the fraction of immune cells between the HF and non-HF groups.Resting NK cells were significantly upregulated in the HF group.*p < 0.05.

Figure 6 :
Figure 6: Profiling of seven IRHGs associated with HF and evaluation of the ability of CXCR5 and FOS to predict long-term HF.(a-g) Comparative analysis of the expression levels of seven IRGs between the HF and non-HF groups.(h and i) ROC analysis of CXCR5 and FOS.*p < 0.05.

Figure 7 :
Figure 7: (a) The macrograph of hearts in the sham and MI group.(b and c) The representative images of cardiac sections from the sham and MI groups (scale bar = 500 μm) subjected to Masson's staining.(d and e) M-mode echo-cardiograms of the sham and MI groups.

Figure 8 :
Figure 8: (a-c) The mRNA expression levels of Fos, Cxcr5, and Nppa in the sham (n = 5) and MI groups (n = 5).(d and e) The ratios of HW/TL and liver weight (LW)/TL in the sham (n = 5) and MI groups (n = 5).(f and g) The Nppa protein levels in the sham (n = 5) and MI groups (n = 5) were evaluated using western blotting.*p < 0.05; **p < 0.01.