Investigating molecular markers linked to acute myocardial infarction and cuproptosis: bioinformatics analysis and validation in the AMI mice model

Data acquisition

A microarray dataset (GSE48060) related to AMI was obtained from the GEO website database (GEO, https://www.ncbi.nlm.nih.gov/geo/) using the “GEOquery” R package (version 4.1) (Davis & Meltzer, 2007). The GSE48060 dataset (GPL570-55999 platform) including 21 healthy and 31 AMI blood samples were selected for further analysis. Robust multiarray average (RMA) was used to normalize the raw gene expression profiles of the GEO dataset (“affy”, version 4.1.0 of the R package).

Animals and procedures

For animal models of MI, male C57BL/6 mice at the age of 6–7 weeks were purchased from Ningbo University. The study was approved by the use committee at Ningbo University School of Medicine (11649) and conducted in accordance with the National Institutes of Health Guide. Animals were kept at 23 ± 2 °C with 50 ± 5% humidity. Mice were randomly split into two groups (n = 10 each): sham and MI. Randomization was done using the Rand() function in Excel based on weight and cage position. Select one cage at random from all the cages in each group. Two animals are taken from each individual recombination and given a permanent number in the cage. The cages were then randomly assigned to the experimental group. The specific steps of MI mouse model establishment are in the Supplemental file S4. The experimental mice were euthanized according to the experimental endpoint set in the experiment. We sacrificed mice after 3 h, 48 h, and 1 week from the beginning of the animal model. Stack experimental animals in the euthanasia box, and introduce carbon dioxide to the animal. There were no mice left in this study and were used for experiments. After 3 h, 48 h, 1 week of MI, the mice were anaesthetized and sacrificed for collection of heart tissues and follow-up examinations.

Echocardiographic analysis

Echocardiography was performed on mice 2 weeks after MI or sham surgery. The mice were anesthetized with isoflurane and their cardiac function was evaluated using a small animal ultrasound system (Vevo 2100; VisualSonics, Toronto, Canada). Cardiac function indicators, such as left ventricular ejection fraction and fractional shortening, were calculated and averaged from at least three cardiac cycles by a single technician in a blinded manner.

Immunofluorescence

Heart sections were deparaffinized, rehydrated, and subjected to antigen retrieval using 10 mM citrate buffer (pH 6.0) at 95 °C for 10 min. After washing with PBS, the sections were permeabilized with 0.5% Triton X-100, blocked with 5% BSA, and incubated with primary antibodies overnight at 4 °C, followed by incubation of secondary antibodies. DAPI was used to label cell nuclei before imaging with a Leica SP8 confocal microscopy.

The primary antibodies used were as follows: ZNF302 Polyclonal Human Antibody, 1:100, (Thermo Fisher Scientific, Waltham, MA, USA); CBLB Monoclonal Mouse Antibody, 1:100, (Abcam, Cambridge, UK); Abclonal, FITC Goat Anti-Rabbit IgG (H+L) (AS011) 1:100; and FITC Goat Anti-Mouse IgG (H+L) (AS001), 1:100.

Quantitative reverse transcription-polymerase chain reaction

The mouse was sacrificed and serum was collected. Total RNA was extracted using the Trizol method. Reverse transcription was performed with a kit from TransGen Biotech (Beijing, China). The cDNA was amplified using real-time PCR and normalized to the housekeeping gene GAPDH. The CT values were analyzed using Applied Biosystems Software, and relative mRNA expression was calculated using 2^−ΔΔC. The sequences of the primers are following: CBLB (5′-3′): Forward GGTCGCATTTTGGGGATTATTGA, Reverse TTTGGCACAGTCTTACCACTTT; ZNF302 (5′-3′): Forward ATCCTCTATTCCTGTTCCTGGGA, Resverse CAGGCCAAATAACATAGGGAGAC.

2,3,5-triphenyltetrazolium chloride stain

Following a myocardial infarction event, the animals underwent anesthesia, intubation, and thoracotomy. The hearts were excised, and the left ventricles were then precisely sliced into 1 mm thick sections. Heart sections were then incubated with a TTC solution (T130066-5 g; Aladdin, Shanghai, China).

Assessing immune cell infiltration

We used a combination of the CIBERSORT algorithm (https://cibersort.stanford.edu/) and Leukocyte signature matrix22 (LM22) to generate this signature matrix and estimate the relative abundance of 22 types of immune cells in each sample. CIBERSORT calculates inverse fold product p-values using Monte Carlo sampling. Samples with p-values less than 0.05 were considered accurate immune cell fractions. On the basis of the sum of the 22 proportions, each sample contained one immune cell (Newman et al., 2015). Reference data sets we put in supplementary data.

Immune cell infiltration and CRG correlation analysis

In order to further demonstrate the association between CRGs and AMI-associated immune cell properties, CRG expression was correlated with immune cell percentage. Statistical significance was determined using the Spearman correlation coefficient when the p-value was less than 0.05. Corrplot (version 0.92) was then used to display the results.

Clustering of AMI patients without supervision

Initially, a total of seven CRGs were obtained according to the previous study by Tsvetkov et al. (2022). The 31 AMI samples were classified into different clusters by unsupervised clustering (Consensus Cluster Plus, version 4.2.1) using seven CRG expression profiles and the k-means algorithm with 1,000 iterations (Langfelder & Horvath, 2008). The maximal subtype number k (k = 2) was chosen, and the cluster number was determined using the cumulative distribution function (CDF) curve, consensus matrix, consistent cluster number, and consistent cluster score (p > 0.9).

Gene set variation analysis

Gene set variation analysis (GSVA) was performed using the R package “GSVA” (version 4.2.1) to examine the differences in enriched gene sets between different clusters of CRGs. The “c2.cp.kegg.symbols.gmt” and “c5.go.symbols.gmt” files were obtained from the MSigDB website database for further GSVA. GSVA scores were compared between different CRG clusters using the “limma” R package (version 4.2.1) and were considered significant if the |t value of GSVA score| was more than 2.

Weighted gene co-expression network analysis

Weighted gene co-expression network analysis (WGCNA) was performed to identify co-expression modules using the R package of “WGCNA” (version 4.2.1). To ensure accuracy, we used the top 25% of genes with the highest variance for subsequent WGCNA analyses. To construct a topological overlap matrix (TOM), we constructed a weighted adjacency matrix with an optimal soft power. As the minimum module size was set to 100, dissimilarity measures (1-TOM) were computed based on the hierarchical clustering tree algorithm. Modules were assigned random colors. In each module, a module eigengene represented global gene expression profiles. Modular significance (MS) demonstrated the relationship between modules and disease. Clinical phenotype and gene significance (GS) have been described as correlated processes in biology.

Machine learning-based prediction model construction

The “caret” R packages (version 6.0.93) were used to establish machine learning models (random forest model (RF), support vector machine model (SVM), generalized linear model (GLM), and Xtreme Gradient Boosting (XGB)) based on two different CRG clusters. As an ensemble machine learning approach, RF utilizes multiple independent decision trees to predict classifications or regressions (Rigatti, 2017). The SVM algorithm generates a hyperplane in the characteristic space that distinguishes between positive and negative instances with a maximum margin (Tang, Tian & Yang, 2018). A GLM is a multivariate linear regression method that can be used to analyze relationships between categorical or continuous factors and normally distributed dependent factors (Gold & Sollich, 2003). Based on gradient boosting, XGB compares classification error and complexity by comparing boosted trees (Sheridan et al., 2020). DEGs were used as explanatory variables in this study in conjunction with distinct clusters as response variables. The 31 AMI samples were randomly classified into a training set (70%, N = 22) and a validation set (30%, N = 9). Grid search was used by the caret package to tune parameters in these models, and all models were evaluated with default parameters using the 5-fold cross-validation method. Using the “pROC” R package (version 1.18.0), we visualized the area under ROC curves. Consequently, the top five predictive genes were determined to be the most important variables for predicting AMI. For verification of the model’s diagnostic value, ROC curves were examined in datasets (GSE123342 and GSE66360).

Model construction and validation for nomograms

The “rms” R package (version 6.2.0) was used to develop a nomogram model to identify AMI clusters. Scores are assigned to each predictor, and the “total score” represents the sum of the scores. Decision curve analysis (DCA) is a method for evaluating clinical predictive models, diagnostic tests, and molecular markers (Ji & Xue, 2020). Nomogram models were evaluated using the calibration curve and DCA.

Drug-therapy analysis

Drug data were downloaded from the Drug Gene Interaction Database (DGIdb) website (https://dgidb.org/search_interactions). We used two hub genes to predict the underlying drug.

Statistical analysis

All data are presented as the mean ± standard deviation. GraphPad Prism 8.0 software (Version X; La Jolla, CA, USA) was used to analyse the data. Unpaired two tailed Student’s t tests and one-way ANOVA were used for data analysis followed by post hoc tests. p < 0.05 was considered statistically significant.

Activation of immune responses in AMI patients due to dysregulation of cuproptosis regulators

To elucidate the biological functions of regulatory factors in the occurrence and development of AMI, the GSE48060 dataset was used to systematically evaluate the differences in expression profiles of 46 CRGs between AMI and non-AMI patients. Figure 1 shows a detailed flow chart of the research process. A total of seven CRGS were identified as differentially expressed cuproptosis genes. Among them, the expression levels of metal regulatory transcription factor 1 (MTF1) and antioxidant 1 copper chaperone (ATOX) were higher, while the expression levels of mitogen-activated protein kinase kinase 1 (MAP2K1), Phosphodiesterase 3B (PDE3B), lipoyltransferase 1 (LIPT1), glutaminase (GLS), and dihydrolipoamide branched chain transacylase E2 (DBT) genes in AMI serum were much lower than those in the non-AMI control group base on log2 fold change (Figs. 2A–2C). Our next step was to examine whether cuproptosis regulatory factors played a fundamental role in the progression of AMI. Surprisingly, some cuproptosis regulators, such as DBT and GLS, showed strong synergies (coefficient = 0.84). At the same time, LIPT1 and MTF1 showed obvious antagonism (coefficient = 0.59). Furthermore, ATOX1 and DBT were associated with other regulators (Fig. 2D) based on the correlation profiles of these CRGs. To clarify whether there were differences in the immune system between AMI and normal controls, immune-infiltration analysis was performed to determine the proportion of 22 immune cell types infiltrated between AMI and non-AMI control subjects according to the CIBERSORT algorithm (Fig. 2E). The results showed that there were significant differences in T cells CD4 memory activated, NK cells resting, macrophages M0, and neutrophils in patients with AMI (Fig. 2F), indicating that the alternation of the immune system may be the main cause of AMI. At the same time, correlation analysis showed that B cells memory, dendritic cells resting, macrophages M0, mast cells activated, monocytes, neutrophils, NK cells resting, T cells CD4 memory resting, T cells CD4 native, and T cells gamma delts were all associated with cuproptosis regulatory genes (Fig. 2G). These results suggested that CRGs may be the critical factors involved in regulating the molecular and immune infiltration status of AMI patients.

Identifying cuproptosis clusters in AMI

To explore the expression patterns associated with cuproptosis in AMI, we used a consensus clustering algorithm to group 31 AMI samples according to the expression profiles of 46 CRGs. A total of seven CRGs were determined as differentially expressed cuproptosis genes. When the k value was set to 2 (k = 2), the cluster number was the most stable, and the consensus index of CDF curve ranged from 0.2 to 0.6 (Figs. 3A, 3B). When k = 2 to 6 (k = 2–6), a difference between the two CDF curves could be seen in the area under the CDF curve (k and k-1) (Fig. 3C). Additionally, only k = 2 resulted in a consistency score of >0.9 for each subtype (Fig. 3D). Combined with the consensus matrix table, 31 AD patients were finally divided into two clusters, Cluster 1 (n = 23) and Cluster 2 (n = 8) (Table S1). According to tSNE analysis (Fig. 3E), there were significant differences between the two clusters.

Characteristics of cuproptosis clusters and cuproptosis regulators

A comprehensive evaluation of the expression differences between Cluster 1 and Cluster 2 was conducted first in order to investigate the molecular characteristics between clusters. Different CRG expression levels were assessed between the two cluster modes. Cuproptosis Cluster 1 showed high expression levels of PDE3B, LIPT1, GLS, and DBT, while cuproptosis Cluster 2 was characterized by high expression levels of MTF1, MAP2K1, and ATOX1 (Figs. 4A, 4B). In addition, immune infiltration analysis showed that Cluster 1 and Cluster 2 had a high proportion of CD8+T cells, T cells CD4 live, T cells CD4 memory resting, T cells CD4 memory activated in cuproptosis (Figs. 4C, 4D). The abundance of monocells and Neutrophils in Cluster 2 was relatively high. This suggests that cuproptosis Cluster 1 may have a more dominant level of immune infiltration. Using WGCNA, we established co-expression networks and modules for normal and AMI subjects to identify hub genes. GSE48060 variance was calculated for each gene expression, and the top 25% of genes with the highest variance were further analyzed. Soft power value 12 and no scale R led to the identification of co-expressed gene module 2 equal to 0.9 (Fig. 5A). The dynamic cutting algorithm was used to obtain five different co-expression modules with different colors, and the heat map of the TOM was constructed (Figs. 5B–5D). The modular-clinical feature (control and AMI) co-expression was subsequently analyzed using these genes in the five color modules. AMI was most strongly related to the blue module, which included 796 genes (Fig. 5E). Further, yellow module genes correlated positively with mode-related genes (Fig. 5F). In addition, we used the WGCNA algorithm to analyze the key gene modules closely related to cuproptosis clusters. We screened β = 17 and R² = 0.9 to construct the scale-free network (Fig. 6A). In particular, 12 modules containing 5,413 genes were identified as important modules, and the heat map illustrated the TOMs of all genes related to these modules (Figs. 6B–6D). Analysis of modular-clinical features (Clusters 1 and 2) showed that the amethyst module (1,115 genes) was associated with AMI clusters (Fig. 6E). As a result of correlation analysis, the gene for the amethyst module was associated with the selected module (Fig. 6F).

Cluster-specific DEG identification and functional annotation

A total of 188 cluster-specific DEGs were identified by analyzing the intersection of cuproptosis module related genes with AMI and non-AMI individual module-related genes (Fig. 7A). GSVA was used to further explore functional differences between the two clusters related to cluster-specific DEGs. The results showed that Mismatch repair, DNA replication, cell cycle, pyrimidine metabolism, RNA polymerase, nucleotide excision repair, and purine metabolism were all present in Cluster 2. However, ECM preceptor interaction and FC epsilon RI signaling pathway were up-regulated in Cluster 1 (Fig. 7B). Therefore, a number of immune responses and metabolic processes may be mediated in Cluster 2.

Machine learning model construction and assessment

We developed four validated machine learning models, RF, SVM, GLM, and XGB, based on the expression profiles of 52 cluster-specific DEGs in the AMI cohort. To interpret the four models and plot the residual distribution for each model in the test set, we applied the “DALEX” package. The GLM machine learning model had relatively low residuals (Figs. 8A, 8B). On the basis of five-fold cross validation, we calculated the ROC curve for four machine learning algorithms. The ROC curve showed the maximum area under general linear model (GLM, AUC = 0.870; RF, AUC = 0.759; SVM, AUC = 0.648; XGB, AUC = 0.741, Fig. 8C). Based on these results, the GLM showed the best ability to discriminate between patient clusters. Finally, the two most important variables in the GLM (ZNF302, CBLB) were selected as predictive genes for further analysis. The GLM model was then evaluated using a column graph to estimate 31 AMI patients (Fig. 9A). Line graph models were evaluated using calibration curves and DCA. AMI cluster risks and predicted risks are very close based on the calibration curve (Fig. 9B). The DCA indicates the high accuracy of our column graph, which may provide a basis for clinical decision-making (Fig. 9C). ROC curves showed satisfactory performance of the 2-gene prediction model with an AUC value of 0.856 in the GSE123342 (Fig. 9D) and 0.719 in the GSE66360 dataset (Fig. 9E), indicating our diagnosis model is equally efficacious in distinguishing AMI from normal individuals. DGIdb was used for drug prediction of the hub genes (Table S2). There was a positive finding for Fasudil and Fluorouracil, but both of these drugs are new and need further research before they can be used clinically.

Expression of the two genes in AMI tissues

To further verify the expression of these two genes in vivo, 6-week-old C57BL/6 mice were divided into sham and AMI groups. Considering the crucial time intervals associated with myocardial infarction, we established three specific time points for our study. The onset of acute myocardial infarction occurs at the 3-h mark, while the 48-h mark represents the timeframe beyond which surgical intervention for acute myocardial infarction is no longer optimal. During this period, we examined the expression of genes related to myocardial infarction. Additionally, we utilized mice with myocardial infarction for a duration of 1 week to investigate changes in molecular expression levels associated with cuproptosis as myocardial infarction progresses. A schematic diagram shows the animal experiment design (Fig. 10A). Additionally, the results of TTC staining were used to verify the degree of ischemic infarction in mice (Figs. 10B, 10C). Based on the statistical data pertaining to myocardial damage, TTC staining revealed a degree of myocardial damage after a duration of 3 h, albeit not as pronounced as the myocardial damage observed after 48 h and 1 week (Fig. 10C). Additionally, echocardiography was employed to identify alterations in cardiac function among mice. The findings indicated that the ejection and contractile functions of the heart exhibited no significant reduction 3 h post AMI. However, a notable decrease was observed at the 48 h and 1 week marks following MI, as illustrated in Figs. 10D–10F. The serum of mice was collected at three different time points for analysis, revealing an increase in the expression of CBLB and ZNF302 at 3 h post-myocardial infarction. The expressions of CBLB and ZNF302 in serum collected at 48 h and 1 week post-myocardial infarction were significantly higher compared to the sham operation group (Figs. 10G, 10H). The immunoconfocal technique was employed to examine the fluorescence staining outcomes of myocardial sections from mice. Our findings indicate that there was no notable disparity in the expression of the two proteins in mice 3 h following MI (Fig. 10I). Nevertheless, disparities in both markers were observed in the myocardial sections of mice 48 h after MI, revealing an increase in marker expression within the tissue subsequent to MI (Fig. 10I). Notably, the disparity between these two expressions was more pronounced in mice with continuous MI for 1 week (Fig. 10I). These findings suggest that the expression difference of these two molecules is not prominent in the early stages of myocardial infarction, but becomes more pronounced over time. These two molecules may serve as potential diagnostic markers for myocardial infarction, particularly left ventricular infarction. Further research is needed to elucidate their role in the treatment of myocardial infarction and the specific mechanism by which they exert their immune effects. These results suggest that ZNF302 and CBLB are risk factors for the development of AMI in mice, although its specific mechanism in AMI patients requires further study.