Microarray data used for this study included the GSE57338 dataset (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE57338), which included 313 total samples (177 diseased samples, 136 control samples). FRGs were identified using established gene sets including PMID: 33867820 (n = 103), PMID: 33767582 (n = 386), PMID: 32760210 (n = 60), and the FerrDb database (n = 369), yielding 449 total FRGs.

The 449 identified FRGs were used to generate a protein-protein interaction network (https://string-db.org), with the MCC algorithm in the CytoHubba v0.1 (New York City, NY, USA) plug-in being used to select hub genes from this network within the Cytoscape application for further downstream analyses.

Functional enrichment analyses for different HF immune-related subtypes were assessed based on enrichment scores for particular gene sets in individual samples as determined with a non-parametric unsupervised gene set variation analysis (GSVA) approach with the R GSVA v1.38 package (Bioconductor, Boston, MA, USA), which enables the transformation of gene expression data from individual genes in a characteristic expression matrix. Differential expression analyses were then performed with the R limma package (Bioconductor, Boston, MA, USA), with a false discovery rate (FDR) $\le$0.05 being indicative of a significant difference.

There are several tools to calculate immune cell infiltration, such as CIBERSORT, XCELL, SSGSEA, TIMER, etc. To avoid the limitations of the database and to improve the accuracy of the data analysis results, we applied here at least three tools to calculate immune cell infiltration, comparing immune cells between subtypes, immune response gene sets, and HLA gene differences.

Hub modules of interest were introduced into the DGIdb database to facilitate the identification of candidate gene-drug interactions, with the ggplot2 v3.3.4 R package (Hadley Wickham, London, UK) being used to generate a Sankey map with geom alluvium as a means of highlighting the associations between hub genes and drug candidates.

Primary antibodies specific for HDAC1, LNPEP, PSMA1 and PSMA6 were from ABclonal (A19571, A11959, A3460, and A2188), while antibodies specific for $\beta$-actin as well as secondary goat anti-rabbit IgG were from Wuhan Servicebio Technology Co., Ltd. (GB15003 and GB23303, Wuhan, China). A test kit for measuring ferrous ion concentration (E-BC-K773-M) was purchased from Elabscience (Wuhan, China), and the kits for measuring malondialdehyde (MDA) and glutathione (GSH) levels were obtained from Nanjing Jiancheng Institute of Biological Engineering (A003-4-1, A006-2-1, Nanjing, China).

C57BL/6 mice (specific pathogen-free 8-week-old males, 20–25 g; Experimental Animal Center of Nantong University) were housed under standard conditions (21 $\pm$ 1 °C, 55–60% humidity) with ad libitum food and water. All animal research was approved by the Animal Ethics Committee of Nantong University in accordance with the Guide for the Care and Use of Laboratory Animals. Animals were randomly assigned to sham or transverse aortic constriction (TAC) surgical treatment groups (n = 6 each), with TAC-induced pressure overload then being modeled as reported previously. Following this surgical procedure, mice were monitored for 6 h, after which they were returned for standard feeding and care to the experimental animal center. In mice that underwent sham procedures, sutures were passed through the soft tissue located underneath the aortic arch without any corresponding ligation, while all other procedure steps were identical to those in mice that underwent the TAC procedure. Mice were monitored daily after surgical treatment and continued to have ad libitum food and water access.

Over 80% of the mice in both surgical groups survived, with comparable rates in both groups. Surviving mice underwent echocardiographic assessment 8 weeks after surgery, followed by further experimentally appropriate examination.

Following isoflurane-mediated anesthetization, mice were fixed in the supine position on an operating table. The skin of the anterior thoracic region was then prepared, coated with coupling agent, and a high-resolution small animal ultrasound imaging system (Vevo 2100, Visual Sonics, Canada) was used to conduct echocardiographic analyses of these animals using the FMS-250 probe (depth: 2.0–2.5 cm, frequency: 13–24 MHz). Under M-mode ultrasonographic visualization of cardiac structures in these animals, parameters measured using the left ventricular (LV) short-axis and parasternal LV long-axis views included LV internal diastolic diameter (LVIDd), LV internal diameter end-systole (LVIDs), and interventricular septal dimension in diastole (IVSd). LV systolic function-related indices were then computed based upon these values over an average of 3–5 cardiac cycles, including ejection fraction (EF) and fractional shortening (FS).

Mouse hearts were collected and fixed in 4% paraformaldehyde for 48 hours and embedded in paraffin. Next, the tissue was cut into 5 mm cross-sectional sections. The sections were stained with hematoxylin-eosin (HE) to observe the histopathological changes in the myocardium, WGA staining to measure the cross-sectional area of cardiomyocytes, and Masson staining to evaluate the extent of interstitial fibrosis in the mouse myocardium.

After collection, cardiac tissue samples were fixed for 10min with 1% glutaraldehyde, followed by further fixation for 1 h at room temperature using 3% glutaraldehyde. Samples were then further fixed for an additional 1 h at room temperature using 1% osmium tetroxide in 0.1 mol/L dimethylarsinate buffer, followed by embedding in resit. Polymerized ultrathin sections were then stained using uranyl acetate and lead nitrate, followed by TEM imaging to assess mitochondrial morphology and transverse muscle arrangement.

Cells were treated as above, and a bicinchoninic acid (BCA) kit was used to detect protein concentrations in sample lysates. A ferrous ion colorimetric test kit (E-BC-K773-M, Elabscience, Wuhan, China) was then used to detect iron ion levels based on provided directions. Briefly, the provided probe bound to ferrous ions released by cells, with the resultant absorbance at 593 nm being quantified to measure intracellular Fe${}^{2+}$ content.

After cells had been harvested and lysed, a BCA kit was used to measure protein concentrations. MDA levels in these cells were assessed using an MDA assay kit, with absorbance at 530 nm ultimately being measured. Similarly, GSH levels in these cells were analyzed based on the directions provided with a micro-reduced GSH assay kit. Absorbance was assessed at 405 nm, and GSH levels were established based on a standard curve.

RNA was isolated from cardiac tissue samples and cardiomyocytes with Trizol (Servicebio, Wuhan, China). HiScript II Q RT SuperMix (Vazyme, Nanjing, China) was used to prepare cDNA, after which ChamQ Universal SYBR qPCR Master Mix (Vazyme, Nanjing, China) and a fluorescent qPCR instrument (Bio-rad; CFX, CA, USA) were used for qPCR analyses. Primers used for this study (RiboBio, Guangzhou, China) had the following sequences: GAPDH: (F: CCTCGTCCCGTAGACAAAATG, R: TGAGGTCAATGAAGGGGTCGT); HDAC1: (F: TCATAAATACGGAGAGTACTTCCCA, R: CGTCTCGCAGTGGGTAGTTCAC); LNPEP: (F: CAAACCAGTCAGCAGAACTCATC, R: TAGAGGAATAATGGCAGTGGGAAG); PSMA1: (F: CGGGTCTAACTGCTGATGCC, R: CATATCGCTGTGTTGGGATCTG); PSMA6: (F: GCAGAGATTGACGCTCACCTT, R: GGCCTCTATGATGTCAGTTTGTTG). GAPDH was used as a control to normalize gene expression, with relative expression being measured using the 2${}^{-\mathrm{\Delta }\mathrm{\Delta }CT}$ method.

A protein extraction reagent was used to isolate proteins from tissue or cell samples, after which a BCA protein analysis kit (Servicebio, Wuhan, China) was used to measure protein concentrations. Proteins were then separated via 8% SDS-PAGE (Servicebio, Wuhan, China), transferred to polyvinylidene difluoride (PVDF) membranes, and blots were then blocked for 2 h at room temperature with 5% non-fat milk prior to incubation overnight with appropriate primary antibodies at 4 °C. Blots were rinsed with TBST three times, then probed with secondary antibodies for 2 h. Chemiluminescent reagents were then used to detect protein bands, followed by quantification using ImageJ V1.8.0 (National Institutes of Health, Bethesda, MD, USA).

Data are reported as means $\pm$ SEM, and were analyzed using GraphPad Prism 9.3.0 (GraphPad Software, San Diego, CA, USA). Differences between groups were compared via Student’s t-tests or one-way ANOVAs with Bonferroni post hoc tests. p 0.05 was the significance threshold.

To gain initial insight into the relationships between ferroptosis and HF, a comprehensive set of 449 FRGs with corresponding genomic location information was established (Fig. 1A). Protein-protein interaction analyses were also conducted for these FRGs (Fig. 1B). The relationship between these FRGs was then examined by comparing all samples in the selected microarray dataset to samples from HF patients (Fig. 1C). Differential FRG expression between healthy control individuals and HF patients was then compared using established cutoff criteria, revealing 62 differentially expressed FRGs of which 42 and 20 were respectively upregulated and downregulated in HF samples (Fig. 1D–F, Supplementary Table 1).

Fig. 1.

Ferroptosis-related genes (FRGs) landscape in heart failure. (A) Genomic position configuration of 449 FRGs. (B) Protein-protein interactions of 449 FRGs. (C) Differential FRG expression correlations (Expression correlation analysis in all samples (top) and in disease samples (bottom), scatter plots showing expression of positively correlated FRGs in different samples). (D) volcano plot showing up-or down-regulation of FRGs in disease compared to control. (E) Box plot showing significant differences in genes in the samples. (F) Heat map showing the expression of FRGs in all samples.

The immunological microenvironment of HF-related patient samples was next characterized. Analyses of the HLA expression profiles in these samples revealed significant downregulation of most HLA genes in HF samples relative to healthy tissues (Fig. 2A). Correlations between these HLA genes and the 62 identified differentially expressed FRGs were also examined (Fig. 2B). Immune infiltration analyses performed through three computational approaches revealed higher levels of B cell and CD8${}^{\mathrm{+}}$ T cell infiltration in HF tissues (Fig. 2C). Comparative analysis of immune-associated genes further indicated that certain cytokine-encoding genes were upregulated in HF samples, including the genes encoding IL-2, IL-6, IL-11, IL-18, IL-21, and IL-22 (Fig. 2D). Correlations between these immune-related genes and HF-related FRGs were also examined (Fig. 2E).

Fig. 2.

The role of FRGs in the immune regulation in heart failure. (A) Expression levels of HLA-related genes in disease and control samples. (B) Correlation analysis of HLA gene expression with FRGs in disease. (C) Immune cell infiltration analysis by XCell, CIBERSORT, and EPIC methods. (D) Expression levels of immune-associated genes in disease and healthy samples. (E) Correlation analysis of immune-associated genes and FRGs in disease samples (ns, not significant, *: 0.05, **: 0.01, ***: 0.001, ****: 0.0001).

A consensus clustering strategy was next used to group HF patient samples based on the 62 identified HF-associated FRGs (Fig. 3A). Samples were separated into two groups based on the fact that the strongest intra-group and weakest inter-group connections were evident at a clustering variable (k) value of 2 in generated CDF curves (Fig. 3B,C). A subsequent principal component analysis (PCA) indicated that these two HF subtypes were distinct from one another based on gene expression profiles (Fig. 3D). Differences in the 62 HF-associated FRGs were also compared between these two subtypes of HF (Fig. 3E,F).

Fig. 3.

Identification of FRGs-related molecular subtypes in heart failure. (A) Heat map showing the results of consistent cluster analysis for subtype identification. (B) The CDF value of consensus index. (C) Relative change in area under CDF curve for k = 2–5. (D) Principal component analysis showing the clustering of each disease sample according to subgroup groupings. (E) Box plot showing the distribution of FRGs expression according to subgroups, with the x-axis representing FRGs and the y-axis the logarithm of relative expression. (F) Heat map showing the expression of FRGs according to subgroups. (ns, not significant, *: 0.05, **: 0.01, ***: 0.001, ****: 0.0001). CDF, cumulative distribution function.

To clarify the biological differences between these HF patient clusters, a GSVA enrichment analysis was conducted for each of these subtypes. KEGG pathway enrichment analyses revealed that relative to samples in Cluster 1, those in Cluster 2 were more closely associated with immune-associated pathways including the primary immunodeficiency, natural killer cell-mediated cytotoxicity, Fc$ϵ$RI signaling, and antigen processing and presentation pathway (Fig. 4A). GO analyses similarly revealed the enrichment of immune-related functions including peptidoglycan immune receptor activity and natural killer cell-mediated immune responses to tumor cells in these HF patient clusters (Fig. 4B).

Fig. 4.

Functional analysis of different molecular subtypes. (A) KEGG. (B) Gene Ontology.

The immunological microenvironment associated with these established HF patient subtypes was next explored by using three algorithms to assess immune cell infiltration in these patient samples. This approach revealed increased immune cell infiltration in samples from Subtype 2 relative to Subtype 1 (Fig. 5A). Moreover, Subtype 2 samples exhibited significant increases in the expression of immune-related genes including IL-6, IL-26, IL-18, TGFB1, and TGFB3 (Fig. 5B). These results highlight key differences in immune cell infiltration among these ferroptosis-associated HF patient subtypes, underscoring the close interplay between ferroptotic processes and the HF-related immune microenvironment.

Fig. 5.

Immune microenvironment of different molecular subtypes. (A) Analysis of immune infiltrating cells between subtypes. (B) Differences in expression of immune response genes between subtypes, ns: not significant, *: 0.05, **: 0.01, ***: 0.001, ****: 0.0001.

A co-expression network was next developed based on the identified differentially expressed genes (DEGs) and HF patient cluster phenotypes in an effort to define gene modules significantly associated with these clusters. Differential gene expression analyses were initially performed among these clusters with established significant thresholds, revealing genes and functions that were distinct among these clusters (Fig. 6A–C). Interactions between these DEGs were then examined through the construction of a co-expression network (Fig. 6D,E). The DynamicTreeCut method revealed three total co-expression modules (Fig. 6F), and the modules that were most highly correlated were filtered based on cluster information. For example, the blue module from Cluster 2 was most strongly correlated with the turquoise module from Cluster 1 (Fig. 6E). Given the large numbers of genes in these modules, they were further assessed to identify hub genes of interest by calculating module feature vectors and their associations with genes to define module membership (MM) and gene significance (GS). In total, 111 hub genes in Cluster1 and 32 hub genes in Cluster2 were identified based on established criteria (GS $>$0.3 and MM $>$0.7) (Fig. 6E,F, Supplementary Tables 1,2).

Fig. 6.

Weighted Coexpression Network Construction and Identification of Hub genes. (A) Volcano plot showing the differential genes between clusters. (B,C) Functional analysis between two clusters. (D) Clustering of samples and removal of outliers. (E) Analysis of network topology for various soft-thresholding powers. (F) Cluster dendrogram of genes in all patients. Each branch in the figure represents one gene, and each color represents one coexpression module. (G) Correlation between the gene module and clusters. (H) The scatter plot of Gene significance vs. Module membership in the turquoise co-expression module. (I) The scatter plot of Gene significance vs. Module membership in the blue co-expression module.

Next, DGIdb was utilized to identify candidate drugs and compounds with the potential to treat HF based on the hub genes identified above. Accordingly, gene-drug interaction information was used to establish a protein-protein interaction network, with hub nodes in this network then being leveraged to determine candidate drugs of interest (Fig. 7A,B). In total, 27 drug candidates associated with 4 hub genes (HDAC1, LNPEP, PSMA1, and PSMA6) were identified (Fig. 7C). Of these compounds, 21 were predicted to interact with HDAC1, including abexinostat, an-9, apicidin, belinostat, cudc-101, and dacinostat.

Fig. 7.

Relationships between core genes and potential drugs. (A) Protein–protein interaction (PPI) network diagram of Cluster1 hub gene. (B) PPI network diagram of Cluster2 hub gene. (C) Sankey diagram showing the potential drugs corresponding to each cluster hub core gene.

To confirm the in vivo relevance of the above findings, an in vivo model of HF was established by performing TAC surgery on C57BL/6 mice, with echocardiographic analyses of these animals being performed at 8 weeks post-surgery (Fig. 8A). These analyses revealed that LV EF and LV FS were significantly reduced in mice in the TAC group relative to the sham control group (Fig. 8B), while the IVSd, LVIDd, and LVIDs were increased (Fig. 8C). Heart volume was also increased significantly in mice in the TAC surgery group (Fig. 8D). HE staining results also showed myocardial hypertrophy in the TAC group of mice, with a large number of disorganized cardiomyocytes, significantly increased cytoplasm, and extensive inflammatory cell infiltration (Fig. 8D). Meanwhile, WGA staining demonstrated a significantly larger ventricular cross-sectional area in TAC mice (Fig. 8E). Subsequently, we assessed myocardial interstitial and perivascular fibrosis in ventricular tissue in Masson-stained sections. As expected, TAC surgery promoted the progression of cardiac fibrosis (Fig. 8F). These results confirmed successful HF model establishment.

Fig. 8.

Validation of the four hub genes expressions in animal models. (A–C) Mice in the indicated groups underwent M-mode echocardiographic imaging to assess left ventricle parameters in the long axis of the left parasternal sternum. Analyzed parameters included ejection fraction (EF), fractional shortening (FS), left ventricular internal diastolic diameter (LVIDd), left ventricular internal diameter end-systole (LVIDs), and interventricular septal dimension in diastole (IVSd). (D) Representative of gross cardiac morphology and H&E-stained cross-sections (20$\times$ and 400$\times$). (E) WGA staining and cell surface area analyses of cardiomyocyte cross-sectional area (CSA). (F) Representative images of interstitial and perivascular myocardial fibrosis in the ventricular sections. (G) Representative images of the ultrastructural morphological characteristics of cardiomyocyte mitochondria. (H) Iron, malondialdehyde (MDA), and glutathione (GSH) levels in murine cardiomyocytes were quantified with appropriate commercial kits. (I) The relative expressions of HDAC1, LNPEP, PSMA1 and PSMA6 were assessed via qPCR in cardiomyocytes. (J) Representative Western immunoblots with corresponding quantification assessing HDAC1, LNPEP, PSMA1 and PSMA6 levels in cardiomyocytes. Data are means $\pm$ SEM. *p 0.05. n = 3/group. Statistical differences between two groups were determined by Student’s t-test. WGA, wheat germ agglutinin; H&E, hematoxylin-eosin; qPCR, quantitative polymerase chain reaction.

Based on the successful establishment of this model, we observed the changes of mitochondria in cardiomyocytes of TAC mice using transmission electron microscopy. The results showed that compared with the Sham group, the mitochondria in the cardiomyocytes of mice after TAC were apparently smaller, disordered in arrangement, with indistinct cristae and ruptured outer membranes (Fig. 8G). Meanwhile, we measured the ferrous ion content and the expression levels of ferroptosis -related factors MDA and GSH, and the results showed that the mice in the TAC group had obviously higher ferrous ion content and significantly increased MDA levels, while GSH levels were significantly lower (Fig. 8H), which further confirmed the close association between HF and ferroptosis.

Subsequent qPCR analysis showed that mRNA expression of HDAC1, LNPEP, PSMA1 and PSMA6 differed significantly between these two groups of mice (Fig. 8I). Western blot analysis also showed differences in the expression of these four genes at the protein level, which was consistent with the results of the bioinformatics analysis described above (Fig. 8J). All of these evidences support a possible relationship between immune microenvironment-related ferroptosis and HF.