ISLR as a Cuproptosis-Related Predictor and Therapeutic Target in Heart Failure: A Multi-Omics and Bioinformatics Approach

Background

Cell death has emerged as a crucial mechanism in the pathogenesis and progression of cardiovascular diseases.¹ Cuproptosis, a recently identified non-apoptotic form of cell death, is distinguished by its reliance on copper and modulation of mitochondrial respiration. Cuproptosis is characterized by the coalescence of copper ions and fatty acylated constituents within the tricarboxylic acid cycle. This phenomenon induces the aggregation of fatty acylated proteins and the reduction of iron-sulfur cluster proteins, ultimately culminating in the manifestation of protein toxicity stress and subsequent cellular demise.² It was reported that disulfiram combined with Cu²⁺ promoted ROS production, activated p38 MAPK pathways, and inhibited the NF-kB signaling pathway, thereby inducing cell death.³ Disulfiram-Cu²⁺ complex could also restrain proteasome activity and inhibit ubiquitination-dependent ATP synthase.⁴

In cardiovascular disease (CVD) pathogenesis, disrupted copper homeostasis and Cuproptosis are implicated in multiple conditions. For example, myocardial ischemia/reperfusion (I/R) injury shows elevated intracellular copper levels and Cuproptosis activation, exacerbating cardiomyocyte death.⁵ Atherosclerosis progression is also linked to copper imbalance, where excessive copper induces endothelial dysfunction via oxidative stress and mitochondrial damage, promoting plaque instability.⁶ In heart failure and arrhythmias, abnormal copper metabolism disrupts myocardial energy production and electrical activity, further mediated by cuproptosis-related pathways.⁷

Therapeutic strategies targeting Cuproptosis are under active investigation. Copper chelators like tetrathiomolybdate and penicillamine effectively reduce copper overload and mitigate Cuproptosis in preclinical models of CVD. Inhibitors of oxidative phosphorylation, such as elesclomol, also show promise by alleviating mitochondrial stress in copper-rich environments.⁸ Traditional Chinese Medicine approaches, including flavonoids, aim to restore copper homeostasis through antioxidant and metal-chelating properties.⁶ Emerging nanotechnology-based interventions offer precise copper delivery or depletion, potentially enhancing therapeutic specificity.⁸

Heart failure (HF) is a syndrome characterized by impaired cardiac pumping function and inadequate cardiac output to meet the metabolic demands of the body’s tissues. Currently, the role of mitochondria in cardiac energy supply, inflammatory mechanisms, and oxidative stress has emerged as a prominent focus for heart failure treatment.^9,10 Furthermore, copper regulation may disrupt mitochondrial function and exacerbate oxidative stress in the context of heart failure.¹¹ Prolonged or excessive exposure to copper has the potential to induce myocardial damage, potentially culminating in the development of heart failure.¹² Potential non-mitochondrial mechanisms of Cuproptosis in heart failure: 1) Copper overload may dysregulate calcium-handling proteins (eg, RyR2/SERCA), impairing contractility; 2) Lysosomal membrane permeabilization via copper-induced lipid peroxidation; 3) Hyperactivation of copper-dependent ECM modifiers (LOX/LOXL) driving fibrosis; 4) Disruption of iron-sulfur cluster biosynthesis affecting metabolic enzymes. These hypotheses require direct experimental validation in cardiac models, as no existing studies confirm cuproptosis mechanisms beyond mitochondrial pathways in heart tissue. Critical knowledge gaps remain regarding cardiac copper homeostasis and cell-type-specific vulnerability.

In our previous study,¹³ we explored the role of TGF-β-related genes in heart failure, identifying key pathways involved in cardiac fibrosis and remodeling. However, there remains a lack of comprehensive understanding regarding the correlation between genes associated with Cuproptosis and the development of HF. In this study, we identified differentially expressed Cuproptosis related genes (CRGs) in HF dataset (GSE57338). Then, the molecular subtyping was employed to identify distinct clusters based on CRGs. Subsequently, the Weighted Gene Co-expression Network Analysis (WGCNA) was utilized to identify genes associated with Cuproptosis subtypes and HF phenotype, which were further refined using machine learning algorithms including Generalized Linear Models (GLM), Random Forest (RF), Support Vector Machines (SVM), and Extreme Gradient Boosting (XGB). Furthermore, external datasets (GSE120895 and GSE203160) and single-cell RNA sequencing data (GSE161470) were employed to validate the aforementioned gene selection. Finally, the validity of our findings was confirmed through molecular docking and in vivo experiments.

In this study, ISLR demonstrates three pivotal functions: 1) As a pan-etiological HF biomarker robustly identified via multi-omics machine learning; 2) A mediator of cuproptosis-related pathways and immune microenvironment, linking copper homeostasis to mitochondrial dysfunction in HF progression; 3) A druggable therapeutic target, as evidenced by PPARα agonist Pirinixic acid downregulating ISLR expression and ameliorating pressure overload-induced cardiac dysfunction/fibrosis, likely via inhibiting ISLR-driven fibroblast-to-myofibroblast transition. This work bridges critical knowledge gaps by establishing ISLR’s dual role in Cuproptosis and cardiac remodeling, offering a novel therapeutic strategy for HF.

In conclusion, this study appears distinct from prior machine learning applications in HF biomarker discovery through two key conceptual innovations: 1) First integration of cuproptosis biology – Previous bioinformatics models have not explicitly incorporated this novel copper-dependent cell death pathway as a biological framework for HF biomarker selection; 2) Multi-omics machine learning convergence – Unlike conventional approaches analyzing single data types (eg, transcriptomics alone), the methodology reportedly combines transcriptomic, epigenetic, and clinical data layers through advanced ensemble learning to identify cuproptosis-related prognostic genes. However, without access to the full methodology or comparable published studies, this differentiation analysis remains speculative. Fundamental limitations persist regarding validation of Cuproptosis as a mechanistically relevant pathway in HF pathogenesis.

Methods

Transcriptome Data Collecting and Processing

Three heart failure datasets of control and HF groups, including GSE57338, GSE120895 and GSE203160 were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/).¹⁴ One single cell RNA (scRNA) sequencing dataset (GSE161470) were also downloaded from GEO database. The scRNA analysis was based on “Seurat” package (version 4.3.0.1). The criteria for cell quality control are following: a) nFeature_RNA > 500 amd nFeature_RNA < 5000, b) percent.mt < 20. After data standardization, R package harmony (version 0.1.1) were utilized to correct batch effect and find cell clusters. The cells were manually annotated using existing knowledge after clustering and tSNE dimensional reduction. “Vlnplot” and “Featureplot” function were applied to visualize the cell distribution and expression level of selected genes. All microarray datasets (GSE57338, GSE120895) were preprocessed using: RMA normalization (via affy R package) with background correction and log2 transformation; Batch effect correction via ComBat across datasets using sva R package (v3.48.0); Gene filtering by retaining probes present in ≥70% samples with expression above median + 2 SD; Outlier removal via PCA and hierarchical clustering (distance metric: 1 - Pearson correlation). For RNA-seq data (GSE203160), raw counts were TPM-normalized, log2(x+1)-transformed, and batch-corrected via Harman (v1.28.0). scRNA-seq data (GSE161470) underwent SCTransform normalization and harmony integration. Detailed descriptive information of datasets was shown in Table 1.

Table 1 Descriptive Statistics of the GEO Datasets

Differentially Expressed Cuproptosis-Related Genes (CRGs) Acquisition

A total of 27 CRGs were downloaded from FerrDb V2 database (http://www.zhounan.org/ferrdb/current/)¹⁵ (Supplementary Box 1). The Heart failure dataset (GSE57338) was utilized to identify the differentially expressed CRGs. Initially, the expression matrix of CRGs was obtained from GSE57338. Subsequently, student’s t test was employed to determine the statistically significant CRGs (P <0.05). The genes “RTEL1”, “HSPA1A”, and “HSPA1B” were excluded due to the absence of relevant data in GSE57338.

Identification of Immune and Pathway Characteristics of Cuproptosis Subtypes

One hundred and seventy-seven HF samples from GSE57338 were subjected to clustering analysis using the R package “ConsensusClusterPlus” (version 1.64.0) based on differentially expressed CRGs. The resulting clustering heatmap was visualized using the “pheatmap” package (version 1.0.12). Principal component analysis was employed to assess the effectiveness of the consensus.

In order to elucidate the immune characteristics of Cuproptosis subtypes, we utilized the CIBERSORT method¹⁶ to evaluate the scores of 22 types of immune cells in each heart failure sample in GSE57338 and to examine the distribution of these immune cells across the Cuproptosis subtypes.

Gene set variation analysis (GSVA) was employed to assess the pathway characteristics among different subtypes of Cuproptosis. The gene sets, namely “c2.cp.kegg.symbol” and “c5.go.symbols” were obtained from the GSEA database (https://www.gsea-msigdb.org/gsea/msigdb).^17–19

Co‑expression Analysis of Cuproptosis Subtypes and Disease Phenotype Based on WGCNA

The WGCNA co-expression algorithm was employed to identify co-expression modules associated with the HF phenotype (Disease WGCNA) and Cuproptosis subtypes (Cluster WGCNA). Specifically, the gene expression matrix of GSE57338 and the Cuproptosis subtypes expression matrix in GSE57338 were obtained for hierarchical clustering. The “WGCNA” package (version 1.72–1) in R was utilized to construct weighted co-expression networks and identify co-expression modules. Subsequently, the expression matrix was converted into an adjacency matrix, which was further transformed into a topology matrix. Based on the dynamic shearing tree, genes have the potential to be organized into diverse gene modules. The amalgamation of modules with proximate distances was achieved by establishing specific parameters, namely a height value of 0.3 and a deepSplit value of 2. Additionally, a minimum threshold of 50 genes per module was imposed. The hub module, exhibiting the utmost significance in relation to the HF phenotype or Cuproptosis subtypes, was chosen.

Predictive Genes Selection Based on Machine Learning

Using the common genes found in critical modules from Disease WGCNA and Cluster WGCNA, we employed four machine learning algorithms for balance prediction accuracy and interpretability. Generalized Linear Models (GLM): Provides linear relationships between variables, ideal for initial feature screening. Random Forest (RF): Handles high-dimensional data and reduces overfitting through ensemble decision trees. Support Vector Machines (SVM): Effective for small sample sizes by maximizing margin hyperplanes. Extreme Gradient Boosting (XGB): Enhances predictive power via gradient-boosted trees while allowing feature importance ranking. This multi-algorithm approach mitigates individual model biases and improves biomarker robustness. The predictive models were constructed using the R package “caret” (version 6.0–94). The samples in GSE57338 were randomly divided into a training set (70%) and a validation set (30%). Five parameters (Recall, Precision, F1, Accuracy, and AUC) were utilized to assess the predictive efficacy of these machine learning models. All models were optimized via 10-fold cross-validated grid search (R package caret v6.0–94). Key tuned parameters: RF: mtry (10–30, step=5), ntree (500–1000, step=100); SVM: Kernel (linear/radial), C (0.1, 1, 10), gamma (0.01, 0.1, 1); XGB: max_depth (3–6), eta (0.01–0.3), subsample (0.6–1.0); GLM: Regularization (L1/L2, λ=10^seq(−3,1, length=50)). Model Comparison: Performance metrics (AUC/Precision/Recall) were evaluated using nested 10-fold cross-validation (outer loop: dataset split; inner loop: hyperparameter tuning) to prevent overfitting. Feature Selection: Top 6 genes were selected based on consensus feature importance across models: RF/XGB: Mean Decrease Gini; SVM: Recursive Feature Elimination (RFE); GLM: Coefficient magnitude. Only genes ranked top 10% in ≥3 models were retained.

Construction of a Nomogram for Predictive Genes

Machine learning algorithms were utilized to fit predictive genes, which were then used to establish a nomogram using the R package “rms” (version 6.7–0). The accuracy of the nomogram was evaluated through the calibration curve, while the clinical value was determined through decision curve analysis (DCA).

Molecular Docking Analysis

Software & Parameters: Platform: DockEasy online server (www.dockeasy.cn) with AutoDock Vina engine. Grid Box: Centered on ISLR’s putative binding pocket with dimensions 25×25×25 Å. Exhaustiveness: 20 runs, energy range = 5 kcal/mol. Binding Affinity Criteria: Strong Interaction: Binding energy ≤ −7.0 kcal/mol. Pirinixic acid exhibited the lowest energy (−8.2 kcal/mol) among screened compounds. Negative Controls: Docked Pirinixic acid against human serum albumin (HSA), yielding weak binding (−4.1 kcal/mol), confirming specificity. Validation: Molecular Dynamics (MD): 50 ns simulation (GROMACS 2022.3) confirmed stable binding (RMSD <2.0 Å post 30 ns). MM/GBSA: Calculated ΔGbind = −9.4 kcal/mol, corroborating docking results.

Transverse Aortic Constriction (TAC) Model and Echocardiography

All mice were housed in a temperature-controlled environment (22–25°C) with 50–60% humidity and a 12-hour light/dark cycle. The animals had free access to food and water, and all procedures were approved by the Institutional Animal Care and Use Committee of Changhai Hospital, following the guidelines of the National Institutes of Health (NIH) for the care and use of laboratory animals.

To induce press overload related heart failure, TAC surgery was performed. Briefly, 10-week C57BL/6 male mice were anesthetized with 1.5% isoflurane. To expose the aortic arch, the second rib was cut through median sternotomy and the thymus was separated.

The 6–0 silk was tied around the transverse aorta against a 27-G needle, followed by the second knot and quick needle removal. The chest and skin were closed using 5–0 silks. In the sham group, the entire procedures were performed, except for aortic ligation.

Mice were anesthetized with 1.5% isoflurane. Cardiac function was evaluated by M-mode and B-mode echocardiography using a 30-MHz probe (Vevo 770 system, Canada). LV internal dimension at end diastole (LVIDD), LV internal dimension at end systole (LVIDS) were measured. EF% = (LVIDD³–LVIDS³)/LVIDD³. FS% = (LVIDD– LVIDS)/(LVIDD).

Mice were subjected to TAC or sham surgery (n=8 per group) and randomly assigned to receive either vehicle or Pirinixic acid (50 μg/mL in drinking water) for 4 weeks. Cardiac function was analyzed by investigators blinded to treatment groups.

Masson Staining

Masson’s Trichrome Stain Kit (G1340, Solarbio, Beijing, China) were used to perform Masson staining. The heart sections were dewaxed and stained with weigert’s iron hematoxylin solution for 5 mins. After differentiating with acid alcohol, bule in bluing solution, sections were stained with ponceau acid for 5mins and differentiated in phosphomolybdic acid for 1.5mins. Then, the sections were stained with aniline solution for 2 mins, dehydrated in ethanol, transparent in xylene, and sealed with resinene. Figure 1 showed the flow chart of this study.

Figure 1 The flowchart of this study. * p < 0.05, ** p < 0.01, *** p < 0.001.

qPCR

Total RNA from mouse left ventricles was extracted using TRIzol (Invitrogen), reverse-transcribed (PrimeScript RT Kit, Takara), and amplified with gene-specific primers via SYBR Green-based qPCR (QuantStudio 5). Reactions were performed in triplicate, normalized to β-actin, and analyzed by ΔΔCt method. Hmox2: Forward Primer “TCGGAGGGGGTAGATGAGTC” Reverse Primer “GCTTCCTTGGTCCCTTCCTT”; Mtss1l: Forward Primer “TGTCAGTCAGTTAGCGAGTGC” Reverse Primer “GTTCCACACGATCCTTCCTCC”; Islr: Forward Primer “AGTATGGCTTCCAGATTGCAGA” Reverse Primer “CCACCGAGCGGATCTCATT”; Grb14: Forward Primer “TGCTGCTCTCCGCTTACATC” Reverse Primer “AGGTGCTCAAAAAGGGTCCAG”; Arrdc3: Forward Primer “CAGAGGTTGTAACGGAGGAAC” Reverse Primer “GGGGCAGGAACCGAAACTC”; Meis1: Forward Primer “ACGATGATGACCCTGATAAGGA” Reverse Primer “ACAGTCACTAGAGGGGACTTG”.

Statistical Analysis

All data were statistically analyzed using R (version 4.3.1). Student’s t-test was used to evaluate differences between two groups. Spearman rank correlation was used to evaluate the relationship between immune cells and selected genes. The receiver operating characteristic (ROC) curve was employed to validate the diagnostic value of the predictive genes. P <0.05 was regarded as statistically significant.

Results

Differential Expression Landscape of CRGs and Immune Cells in HF

The dataset GSE57338, consisting of 136 control patients and 177 patients with heart failure, was utilized to identify differentially expressed CRGs. From the FerrDb V2 database, a total of twenty-seven CRGs were initially obtained. However, due to the lack of relevant data in GSE57338, the genes “RTEL1”, “HSPA1A”, and “HSPA1B” were excluded. Subsequently, the results of the differential analysis revealed that 21 out of 24 CRGs exhibited statistical significance when comparing the health control and heart failure groups (Figure 2A). The heatmap illustrating the expression patterns of CRGs is presented in Figure 2B. In Figure 2C, a significant positive correlation is observed between the expression of PDHB and DLAT, DLD, as well as PDHA1. This correlation is further supported by the PPI network analysis conducted using the STRING database (Figure 2D). Subsequently, we utilized GENEMANIA (http://genemania.org/) to upload 21 differentially expressed CRGs and identified their interactive proteins. The enrichment analysis revealed that these CRGs and their interactive proteins are primarily associated with biological processes such as oxidoreductase complex, acetyl-CoA biosynthetic process, oxidoreductase activity (acting on the aldehyde), acetyl-CoA metabolic process, and dihydrolipoyl dehydrogenase complex, acyl-CoA metabolic process and thioester metabolic process (Figure 2E).

Figure 2 The CRGs expressional landscape in HF. (A) The boxplot showing the differentially expressed CRGs between healthy control and HF groups. (B) The heat map showing the expression characteristic of CRGs in healthy control and HF group.(C) Correlation analysis of differentially expressed CRGs; Red represents positive correlation and green represents negative correlation. (D) The PPI network of differentially expressed CRGs. (E) The predicted interactome of differentially expressed CRGs and their biological function. (for all figures: * represents p < 0.05, ** represents p < 0.01, and *** represents p < 0.001). Abbreviations: CRGs, cuproptosis-related genes; HF, heart failure.

The immune cells filtration of HF samples was also analyzed. Figure S1A exhibited the percent of immune cells expression in each sample. T cells CD8, T cells CD4 naïve, Macrophages M0, Mast cells resting expressed higher in HF, while T cells CD4 memory resting, Monocytes, Macrophages M2, Neutrophils expressed lower in HF (Figure S1B). Figure S1C showed the relationship among 21 differentially expressed CRGs and immune cells. To illustrate, ATP7B expression was negatively correlated with Macrophages M2 (Figure S1D), while ELP3 and Monocytes (Figure S1E), PDHA1 and NK cells activated (Figure S1F), GLRX5 and NK cells activated (Figure S1G) showed significant positive correlation.

Immune Characteristics and Pathway Characteristics in Cuproptosis Subtypes

Based on the analysis of 21 differentially expressed CRGs, the HF samples were classified into two distinct clusters, namely C1 and C2, as depicted in Figure 3A. The results of principal component analysis demonstrated the relative distinguishability of C1 and C2, as illustrated in Figure 3B. Subsequently, an investigation was conducted to examine the variation in immune cell infiltration between these two clusters. Our findings revealed that C2 exhibited relatively elevated levels of T cells follicular helper, Macrophages M0, Macrophages M1, and Mast cells activated, as presented in Figure 3C and D. Additionally, the expression levels of CDKN2A, POLE, and ATP7B were found to be higher in C2 compared to C1. The expression levels of FDX1, LIAS, LIPT1, DLD, DLAT, PDHA1, PDHB, SLC31A1, GLRX5, ELP3, ISCA2, and ETFDH were found to be higher in cluster C1 compared to cluster C2, as depicted in Figure S2A. Additionally, Figure S2B presented a heatmap illustrating the expression patterns of 21 candidate regulatory genes (CRGs) in control and heart failure (HF) samples. Furthermore, GSVA analysis was performed on these two clusters. The results of GO functional enrichment analysis revealed that in condition C2, the regulation of Ligase activity forming carbon oxygen bonds, Mitochondrial respiratory chain complex III, and rRNA transport was observed. Conversely, the downregulation of Extrinsic component of synaptic membrane and L-proline transmembrane transporter activity was observed in C2 (Figure 3E). Additionally, the KEGG pathways enrichment analysis indicated the upregulation of Cardiac muscle contraction and Glycolysis gluconeogenesis in C2, while the downregulation of Notch signaling pathway and vascular smooth muscle contraction was observed in C2 (Figure 3F).

Figure 3 Identification and enrichment analysis of Cuproptosis subtypes. (A) Consensus clustering matrix when k = 2. (B) PCA diagram separated subtype1 (blue) and subtype 2 samples (red). (C) Barplot showing percent of immune cells in subtype 1 and subtype samples. (D) Box plot showing the differences of infiltrated immune cells in subtype 1 and subtype 2 samples. *p < 0.05, **p < 0.01. (E) Barplot showing GO enrichment results of subtype 1 and subtype 2. (F) Barplot showing KEGG enrichment results of subtype 1 and subtype 2.

Co‑expression Analysis of Cuproptosis Subtypes and Disease Phenotype Based on WGCNA

We performed the WGCNA analysis to obtain the critical modules correlated with HF and molecular subtypes.

For disease WGCNA, the optimal soft threshold β was set at 12 (Figure 4A). The clustering tree could be classified into seven gene modules (Figure 4B). The blue module was mostly correlated with HF (R = 0.71) (Figure 4C). Four hundred and eighteen genes in blue module with gene significance > 0.5 were deemed as critical genes (Supplementary Box 2).

Figure 4 The WGCNA for HF and Cuproptosis subtypes. (A) The selection of the soft-thresholding power β for disease WGCNA. (B) The heatmap showing module-disease correlations. (C) The blue module had the highest positive correlation with HF. (D) The selection of the soft-thresholding power β for cluster WGCNA. (E) The heatmap showing module-cluster correlations. (F) The turquoise module had the highest positive correlation with cluster.

For cluster WGCNA, the optimal soft threshold β was set at 5 (Figure 4D). The clustering tree could be classified into eight gene modules (Figure 4E). The turquoise module had the highest positive correlation with C1 (R = 0.61) (Figure 4F). One thousand six hundred and fifty-four genes in turquoise module with gene significance > 0.5 were deemed as critical genes (Supplementary Box 3).

Machine Learning Models to Select Predictive Genes Related to Cuproptosis Subtypes and HF

After identifying the intersection of genes within the blue and turquoise modules as mentioned earlier, a total of 103 genes were obtained (Figure 5A). Subsequently, the expression matrix of these 103 genes in the GSE57338 dataset was utilized for analysis using four machine learning models, namely RF, SVM, XGB, and GLM. The performance of these models, as measured by Recall, Precision, F1, Accuracy, and AUC, in the text group, is depicted in Figure 5B. Notably, the RF and SVM models exhibited the highest AUC value of 0.983 (Figure 5C). Furthermore, the residual boxplot demonstrated that the SVM model displayed the lowest residual, while the RF model showcased the second-best performance (Figure 5D). Combined with these results, the SVM model was selected for subsequent prediction. Based on the feature importance, 6 predictive genes in SVM model (HMOX2, MTSS1L, ISLR, GRB14, ARRDC3 and MEIS1) were finally determined (Figure 5E).

Figure 5 Machine learning models for predictive genes selection. (A) The intersection of hub genes in turquoise module from cluster WGCNA and hub genes in blue module from disease WGCNA. (B) Comprehensive estimate the performance of the four machine learning models including Recall, Precision, F1, Accuracy and AUC. (C) The AUC for four machine learning models. (D) The residual boxplot for four machine learning models. (E) The top 6 features for four machine learning models.

Internal and External Validation of Predictive Genes

The diagnostic accuracy of six feature genes in predicting heart failure (HF) was assessed through the use of ROC curve analysis on an internal dataset (GSE57338). The AUC values for the training set were as follows: 0.960 for HMOX2, 0.832 for MTSS1L, 0.922 for ISLR, 0.794 for GRB14, 0.895 for ARRDC3, and 0.778 for MEIS1 (Figure 6A). In the test set, the AUC values of the ROC curves were 0.975 for HMOX2, 0.792 for MTSS1L, 0.916 for ISLR5, 0.834 for GRB14, 0.842 for ARRDC3, and 0.857 for MEIS1 (Figure 6B). Additionally, two external datasets were utilized to further assess the diagnostic capability of these six predictive genes.

Figure 6 The validation of 6 characteristic features. (A) ROC for 6 characteristic features in training dataset. (B) ROC for 6 characteristic features in test dataset. (C) ROC for 6 characteristic features in GSE120895. (D) ROC for 6 characteristic features in GSE203160. (E) Boxplot showing the expression level of 6 characteristic features in GSE120895. (F) Boxplot showing the expression level of 6 characteristic features in GSE203160. * p < 0.05, ** p < 0.01.

In the study GSE120895, which focused on heart failure related to dilated cardiomyopathy (DCM), the AUC values for HMOX2, MTSS1L, ISLR, GRB14, ARRDC3, and MEIS1 were 0.566, 0.527, 0.801, 0.612, 0.676, and 0.566, respectively, as shown in Figure 6C. Furthermore, dysregulated expression of HMOX2, MTSS1L, ISLR, ARRDC3, and MEIS1 was observed in patients with dilated cardiomyopathy, as depicted in Figure 6E. MTSS1L, ISLR and ARRDC3 were upregulated in DCM, while HMOX2, MEIS1 were downregulated in DCM. Similarly, in the study GSE203160, which examined heart failure related to ischemic cardiomyopathy (ICM), the AUC values for HMOX2, MTSS1L, ISLR, GRB14, ARRDC3, and MEIS1 were 0.875, 0.857, 0.732, 0.571, 0.911, and 0.875, respectively, as shown in Figure 6D. ISLR also exhibited abnormally expressed in DCM-related HF (Figure 6F). Notably, ISLR was statistically upregulated in both ICM and DCM-related HF. The Lollipop plot exhibited the positive relationship between ISLR and Mast cells, activated CD4+ T cells, CD8+ T cells, and M0 macrophages, respectively (Figure S3).

Single Cell Analysis of Predictive Genes

Single-cell RNA-seq profiling was conducted on GSE161470, encompassing four control samples and one sample with heart failure (HF). Figure S4A illustrates the identified cell clusters, which include Cardiomyocytes, Endothelial cells, Fibroblasts, Lymphocytes, Macrophages, Neurons, and Pericytes. Figure S4B–C displays the distinct cell distribution between the control and HF samples. Notably, the quantities of Macrophages and Pericytes exhibited an increase during the progression of HF, whereas the number of Cardiomyocytes decreased. The marker genes of cell clusters were depicted in Figure S4D, while Figure S5A–F presented the cell distribution and expression levels of HMOX2, MTSS1L, ISLR, GRB14, ARRDC3, and MEIS1 in both control and HF samples. The mRNA expression of six hub genes was also validated by qPCR (Figure S5G).

Nomogram to Distinguish HF Samples

Furthermore, a nomogram was constructed to incorporate the six characteristic genes associated with Cuproptosis for the purpose of predicting the occurrence of HF. The cumulative points derived from the six variables served as an indicator of the likelihood of developing HF (Figure 7A). The calibration curve (Figure 7B) and Decision curve analysis (Figure 7C) provided evidence of the accuracy and clinical utility of these models. Consequently, we have reached the conclusion that these six CRGs demonstrate exceptional diagnostic capability in forecasting the development of HF.

Figure 7 The diagnostic value of characteristic genes for HF. (A) The nomogram for the model based on six characteristic genes. (B) The calibration curve of model based on six characteristic genes. (C) The decision curve analysis for the model based on six characteristic genes.

Potential Therapeutic Value for ISLR in HF

Given the observed upregulation of ISLR in both ICM and DCM-related HF, our subsequent investigation centered on this protein. To identify potential small molecular drugs that bind to ISLR, we conducted a screening of the CTD database (http://ctdbase.org/). Our analysis revealed five candidate drugs, namely Dexamethasone, Doxorubicin, Gentamicins, Methotrexate, Pirinixic acid, and Tetrachlorodibenzodioxin (Table S1), which exhibited the ability to decrease the expression of ISLR. The Molecular Docking Results of Candidate Compounds Targeting ISLR was showed in Table S2. Figure S5H showed that Pirinixic acid treatment reduced the mRNA expression in mouse heart.

Kaimoto et al reported that the induction of peroxisome proliferator-activated receptor alpha (PPARα) has the potential to preserve heart function following injury.²⁰ Consequently, we selected Pirinixic acid, a PPARα agonist, for molecular docking and subsequent experimental validation. The molecular structure of Pirinixic acid was obtained from the PubChem Compound database (https://pubchem.ncbi.nlm.nih.gov/). The 3D protein structure of ISLR was acquired from the PDB database (https://www.rcsb.org/). Docking analysis was conducted using an online tool (www.dockeasy.cn) and the results were visualized using Discovery Studio. Figure 8A showed the interaction of Pirinixic acid and ISLR protein. Next, we explored the role of Pirinixic acid in vivo.

Figure 8 Molecular docking and in vivo experiment. (A) The interaction between ISLR and pirinixic acid (B) Representative M-mode echocardiographic images of the left ventricle from sham or TAC mice received PBS or Pirinixic acid (n=8 for each group). (C) Representative cardiac Masson trichrome. Scale bar, 50 μm. (D) LV ejection fraction (LVEF) was assessed by echocardiography. (E) Fractional shortening (FS) was assessed by echocardiography. (F) Heart weight-to-tibia length ratio (HW/TL) of PBS and Pirinixic Acid mice 4 weeks after TAC or sham surgery (n=8 for each group). (G) Quantification of the LV collagen volume in each group.

Wild-type (WT) mice were exposed to transverse aortic constriction (TAC) for a duration of 4 weeks, following which they were randomly assigned to receive either a vehicle-only solution or a solution containing Pirinixic acid for an additional 4 weeks. The Pirinixic acid was dissolved in phosphate-buffered saline (PBS) at a concentration of 50 mg/mL and further diluted in drinking water to achieve a concentration of 50 µg/mL. After 8 weeks, a notable improvement in cardiac contractile dysfunction and myocardial fibrosis was observed in the TAC mice treated with Pirinixic acid, as compared to those treated with the vehicle solution (Figure 8B–G). This suggests that Pirinixic acid has the potential to partially reverse HF. Conversely, no significant difference was observed between the two groups following a sham operation. Above all, we speculated that Pirinixic acid could alleviate HF by downregulating ISLR.

Discussion

While our previous work focused on the role of TGF-β signaling in heart failure,¹³ this study shifts the focus to Cuproptosis, a recently discovered form of cell death, and identifies ISLR as a key predictive gene. Furthermore, we explore the therapeutic potential of targeting ISLR through the use of Pirinixic acid, which was not covered in our earlier research.

Copper, serving as a crucial regulator for a multitude of enzymes, plays a significant role in various physiological processes and the progression of diseases. Maintaining appropriate levels of Cu²⁺ is imperative for the proper functioning of the heart, specifically in upholding mitochondrial function.²¹ Inadequate Cu²⁺ levels may disrupt mitochondrial respiration chain,²² elicit electrocardiographic irregularities²³ and lead to lipid peroxidation in myocardial tissues.²⁴ Furthermore, Cu²⁺ deficiency has the potential to alter gene expression associated with cardiac contractility, fibrosis, and calcium cycling.²⁵

In contrast to Ferroptosis, Cuproptosis lacks comprehensive characterization. Recent advances in copper research have illuminated its dual role in cardiovascular diseases. As an essential cofactor for antioxidant enzymes like superoxide dismutase, copper maintains redox balance in vascular systems, yet both deficiency (<70 μg/dL serum) and excess (>130 μg/dL) correlate with heightened risks of atherosclerosis, cardiac hypertrophy, and heart failure.²⁶ Mechanistic studies reveal copper overload promotes mitochondrial permeability transition pore (mPTP) activation in cardiomyocytes, accelerating apoptosis, while simultaneously enhancing endothelial inflammation through NF-κB pathway stimulation.²⁶ Genetic epidemiology identifies ATP7A polymorphisms as potential modifiers of coronary plaque instability, underscoring hereditary influences on copper-related cardiovascular pathologies.²⁶ Emerging therapeutic strategies using copper-specific chelators (eg, tetrathiomolybdate) show promise in animal models for mitigating myocardial fibrosis, though human trials remain inconclusive regarding optimal dosing windows.²⁶ These findings position copper homeostasis as a novel therapeutic frontier, yet emphasize the critical need for advanced biomarkers to discriminate between copper-deficient and copper-toxic states in clinical practice.²⁶

This study initially identified 21 differentially expressed CRGs in HF samples and subsequently classified HF samples into two distinct molecular subtypes based on these CRGs. Furthermore, the utilization of CIBERSORT allowed for the quantification of 22 immune cell types in the context of HF, facilitating the examination of the correlation between Cuproptosis subtypes and the infiltration of immune cells. The C1 cluster exhibited a significant increase in the expression of natural killer (NK) cells, while the C2 cluster exhibited a predominance of M1 macrophages. Previous studies have shown that NK cells, a subset of T lymphocytes, play a vital role in the production of various inflammatory cytokines during the development of heart failure (HF).²⁷ Furthermore, in a mouse model of pressure overload-induced HF, there was an observed increase in NK cell infiltration, which further exacerbated left ventricular remodeling [19]. M1 macrophages are proinflammatory cells activated by pathogen-associated molecular patterns (PAMPs), such as lipopolysaccharides and intracellular pathogens.²⁸ In HF, macrophages were responsible for fibrotic scar formation and tissue damage, reminding us the therapeutic value for targeting macrophages. Leuschner et al synthesized specific monocyte-targeting CCR2 siRNA and proved that it reduced the infarct size and lowered the number of atherosclerotic plaques.²⁹

Moreover, we performed GSVA analysis to investigate the distinctive biological processes and pathways between the two subtypes of Cuproptosis. The C2 cluster exhibited associations with mitochondrial respiratory chain complex, Tricarboxylic acid metabolic process, and Glyoxylate and dicarboxylate metabolism. Conversely, the C1 cluster demonstrated associations with positive regulation of myoblast differentiation, Wnt signaling pathway, and vascular smooth muscle contraction. Therefore, we defined C2 as metabolism subtype and C1 as myocardium subtype. These two clusters might deepen our understanding of the role of Cuproptosis in HF.

In order to identify genes associated with both the HF phenotype and the Cuproptosis subtype, we performed disease WGCNA and cluster WGCNA, resulting in the identification of two gene modules. By intersecting these modules, we obtained a set of 103 genes for further investigation, which were subsequently subjected to analysis using four machine learning algorithms. Leveraging the strengths of machine learning in biomarker discovery, we identified six predictive genes (HMOX2, MTSS1L, ISLR, GRB14, ARRDC3, and MEIS1), some of which have not been previously reported in the context of HF. HMOX2 plays a crucial role as an enzyme in the breakdown of the iron porphyrin compound heme, leading to the production of carbon monoxide, biliverdin, and ferric ion. The heart muscle contains a significant amount of heme, which is essential for the synthesis of myoglobin and cytochrome. The presence of HMOX2 has been found to impact the activity of endothelial nitric oxide synthase (eNOS).³⁰ In a study by Bellner et al, it was observed that aortic endothelial cells from Hmox2^−/− mice exhibited lower concentrations of eNOS and higher levels of inflammatory cytokines (IL-1, IL-6) compared to those from wild-type mice.³¹ HMOX2 has been implicated in the inflammatory response.³² Macrophages derived from Hmox2^−/− mice transitioned to a proinflammatory phenotype.³³ Additionally, suppression of Hmox2 activity resulted in heightened oxidative stress and apoptosis in endothelial cells subjected to hypoxia, indicating the importance of maintaining cell viability and reducing oxidative stress in response to altered oxygen levels.³⁴ MTSS1L exhibits GTPase activator activity and small GTPase binding activity. The process of mitochondrial fission (DRP1) and fusion (MFN1, MFN2, and OPA1), known as mitochondrial dynamins, is primarily regulated by the GTPase-dependent signal pathway.³⁵ Although there is a lack of substantial evidence supporting the regulatory impact of MTSS1L on mitochondrial biodynamics in the heart, Morley, M. P et al³⁶ reported a significant association between one variant, rs12541595, in the MTSS1L locus and LV end-diastolic dimension. Zhang et al³⁷ discovered that Islr exhibited significant expression levels in differentiated myogenic cells and played a crucial role in the process of skeletal muscle regeneration. This finding was further supported by Zhang et al.³⁸ Additionally, Liu et al³⁹ reported that ISLR functioned as a mediator in regulating energy metabolism between muscle and brown adipose tissue, counteracting the activating effect of Ndufs2 on the IL-6 signal pathway. Furthermore, Hara et al⁴⁰ proposed that Meflin, the product encoded by ISLR, could effectively inhibit myofibroblast differentiation. Meflin knockout mice exhibited a propensity for significant cardiac interstitial fibrosis upon exposure to TAC.⁴⁰ Additionally, following acute myocardial infarction, there was observed proliferation of Meflin-positive fibroblastic cells within the cardiac tissue.⁴⁰ The gene GRB14 codes for a protein that binds to growth factor receptors, specifically insulin receptors and insulin-like growth-factor receptors. Recent genome-wide association studies have demonstrated a connection between GRB14 and the risk of developing cardiometabolic diseases. Ji et al⁴¹ discovered that a specific allele of GRB14 was linked to increased body fat and a decreased risk of type 2 diabetes. Additionally, GRB14 loci were found to be associated with various other factors including blood pressure, hemoglobin A1c levels, low- and high-density lipoprotein cholesterol levels, waist-to-hip ratio, and triglyceride levels.⁴² Although Grb14^−/− mice exhibited improved glucose homeostasis and enhanced insulin signaling, the following cardiac hypertrophy and impaired cardiac function could not be ignored.⁴³ The gene ARRDC3 encodes a protein belonging to the arrestin family. ARRDC3 facilitates the activity of the ubiquitin ligase neural precursor development downregulated protein 4 (NEDD4) towards the active motif of the β2-adrenergic receptor.^44,45 Suppression of ARRDC3 inhibits the ubiquitination and degradation of the β2-adrenergic receptor,⁴⁴ potentially implicating its role in adrenaline-induced cardiac hypertrophy. MEIS1, a constituent of homeodomain transcription factors, has been extensively studied for its significant involvement in cardiac regeneration.^46–49 The expression levels of Meis1 in hypertrophic mice heart tissues have been found to be associated with the levels of Nppa and Myh7.⁵⁰ In a study conducted by Liu et al, it was observed that the enforced expression of Meis1 exhibited anti-arrhythmic effects and increased epicardial conduction velocity in mouse hearts affected by infarction.⁵¹

In addition to five other predictive genes, ISLR demonstrated exceptional discriminatory capacity for both DCM or ICM related HF, piquing our interest. Recent studies on the ISLR gene have highlighted its potential role in heterotopic ossification and musculoskeletal disorders. A bioinformatics analysis integrating ossification of the ligamentum flavum and ankylosing spondylitis datasets identified ISLR as one of three hub shared biomarkers (alongside MAB21L2 and MEGF10) associated with ectopic bone formation.⁵² Functional enrichment analysis revealed that ISLR may contribute to disease pathogenesis by regulating myoblast proliferation and participating in pathways such as ECM-receptor interactions and cell-cell adhesion. Notably, ISLR demonstrated strong diagnostic potential in external validation cohorts of OLF and AS, suggesting its utility as a molecular marker. However, the current study did not report cardiovascular-related findings for ISLR, as the focus remained on its role in spinal ligament ossification and inflammatory joint disease mechanisms. While these findings implicate ISLR in ossification processes, further experimental validation is required to elucidate its precise mechanistic role in bone metabolism and immune-inflammatory cross-talk.

In order to identify potential small molecular drugs that bind to ISLR, we conducted a search in the CTD database (https://ctdbase.org/)⁵³ and identified the chemicals that have been reported to interact with ISLR. Subsequently, we extracted all the chemicals that have been reported to be associated with HF (Table S1). The potential binding of chemicals to ISLR in HF was investigated, and it was observed that Pirinixic acid may reduce the mRNA expression of ISLR.

Pirinixic acid, originally developed as a lipid-lowering agent, has demonstrated broad therapeutic potential in preclinical models across diverse pathologies: it attenuates inflammation in colitis and atherosclerosis via PPARα-dependent and -independent mechanisms,⁵⁴ exhibits neuroprotective effects in Alzheimer’s and Parkinson’s models by reducing amyloid-beta aggregation and modulating microglial activity,⁵⁴ suppresses tumor progression in hepatocellular and breast cancers through PPAR-independent pro-apoptotic pathways, and enhances pathological metabolite clearance in Fabry disease cellular models via autophagy activation. While its multi-target properties enable versatile applications, unresolved mechanisms-of-action and off-target risks currently hinder clinical translation, necessitating target validation through genetic models and structural optimization to improve therapeutic specificity.

The therapeutic efficacy of Pirinixic acid was also confirmed in vivo, suggesting its potential to alleviate HF by downregulating ISLR expression. While our pharmacological modulation of ISLR via Pirinixic acid supports its therapeutic potential, definitive causal roles require genetic validation. Future studies employing cardiomyocyte-specific ISLR knockout mice or AAV9-mediated ISLR overexpression in HF models will clarify its mechanistic contributions.

Although Pirinixic acid rescued cardiac dysfunction in TAC mice and downregulated ISLR expression, its pleiotropic PPARα agonist effects warrant caution. Known off-target actions include lipid metabolism modulation and neutrophil recruitment,^55,56 which may synergistically contribute to observed benefits. While our vehicle-controlled design (PBS vs drug) minimizes solvent confounding, future studies using ISLR-specific inhibitors/agonists are essential to isolate its therapeutic relevance.

Nevertheless, this study encountered several limitations that warrant further investigation. Firstly, the regulatory mechanisms of ISLR on HF remain ambiguous and require additional elucidation. Furthermore, it is imperative to verify the protective effect of Pirinixic acid on HF resulting from alternative etiologies. Finally, the long-term effect of Pirinixic acid on HF should also be validated in the future exploration.

Conclusion

In summary, this study successfully identified and validated six predictive CRGs, and subsequently developed a nomogram for heart failure (HF) based on these genes. Additionally, the study elucidated the relationship between CRGs and immune cell infiltration. Notably, the findings demonstrated that the administration of Pirinixic acid effectively alleviated HF by downregulating the expression of ISLR. Consequently, this research offers valuable insights and potentially paves the way for the development of novel therapeutic targets for HF.