Creation of a Rat Takotsubo Syndrome Model and Utilization of Machine Learning Algorithms for Screening Diagnostic Biomarkers

Introduction

TTS, also known as “stress cardiomyopathy” or “cardiovascular neurosis”, was first discovered by Japanese researchers. This condition is comparable to acute coronary syndrome and is induced by sudden physical or emotional stress. It is often characterized by brief myocardial edema or systolic myocardial dysfunction, which typically takes the form of left ventricular systolic dyskinesia. Acute chest pain, dyspnea, electrocardiographic changes like ST-segment elevation, T-wave inversion, slightly elevated cardiac enzyme levels, and transient ventricular wall motion abnormalities are typical symptoms of acute coronary syndrome.¹ Early imaging is frequently used to distinguish between TTS and acute coronary syndrome,² although accuracy needs to be improved.

Cell death is a crucial physiological process involved in growth, aging, and tissue homeostasis, which is often dysregulated in various clinical situations.³ Ferroptosis, an iron-dependent non-apoptotic form of cell death, was discovered in 2012.⁴ Its mechanism is distinct from those of apoptosis, necrosis, pyroptosis, and autophagy, and it is centered on the accumulation of lipid peroxides, which essentially results from intracellular reactive oxygen species (ROS) metabolism disorders.⁵ Iron catalyzes metabolic irregularities, reducing the antioxidant capacity of the cells, and consequently resulting in a gradual accumulation of ROS that causes cell death.⁶ Ferroptosis also contributes to significant morphological changes in mitochondria, including a reduction in mitochondrial volume, an increase in membrane density, and a disturbance of the mitochondrial membrane’s organization.^7,8 The development of ferroptosis is significantly influenced by circulating iron, which can be inhibited using iron chelators. Furthermore, cells with high transferrin expression on their membranes are more susceptible to ferroptosis.⁴ Iron ions and lipid peroxide levels first increase significantly has been detected in ferroptosis. These events lead to amplified lipid peroxidation reactions on the cell membrane and increase the free radical studies of lipid peroxides, encouraging the apoptotic process and ultimately leading to cardiomyocyte death.⁹ Secondly, numerous studies have shown that TTS may be associated with alterations in the lipid metabolism enzyme system in vivo.¹⁰ The GLS (lipid peroxide-scavenging enzyme) mechanism, which is essential to ferroptosis and may disrupt the balance of lipid metabolism in vivo, is what causes TTS to manifest.¹¹ The degree of glutathione peroxidase (GPX4) degradation and Nrf2 activation in cardiomyocytes, which are associated with the mechanism of ferroptosis, has been observed closely linked to the development of TTS.¹² Taken together, these findings indicates that ferroptosis may play a role in the process of TTS. Mounting evidence suggests that ferroptosis significantly contributes to the development of cardiovascular diseases such as myocardial ischemia-reperfusion injury, hypertension, atherosclerosis, myocardial infarction, and heart failure.^4,13,14 However, the mechanism of ferroptosis in TTS is still unclear. With the advent of Next Generation Sequencing (NGS) technology, bioinformatics has become a burgeoning interdisciplinary field that mainly utilizes molecular biology, computer science, and information technology to elucidate the laws and mysteries underlying complex biological data. Among these, machine learning algorithms have demonstrated a significant role in the screening of sensitive diagnostic biomarkers for various diseases.¹⁵ Intending to provide clinical guidelines for the diagnosis and treatment of TTS, this study conducted a thorough analysis of TTS transcriptome data using machine learning techniques to develop a diagnostic model based on FRGs.

Materials and Methods

Animals

All animal studies were approved by the Southwest Medical University Institutional Review Board’s Animal Welfare Committee and adhered to the institution’s guidelines for animal experimentation (approval number: 20220815-034). Female SPF-grade SD rats and ISO (isoprenaline) were purchased from Beijing HFK Bioscinece Co., Ltd and MedChemExpress respectively for this study. After 24 hours, of intraperitoneal injection of ISO, small animal ultrasound was used to detect the heart function to confirm that TTS rat model was successfully established. Then the rats were euthanized with isoflurane, and heart tissues were collected for RNA-Seq analysis.

RNA Extraction and Library Construction

A total RNA was used to prepare the RNA sample. mRNA was isolated from total RNA using magnetic beads with poly-T oligonucleotides attached. Divalent cations were used to break the RNA at a high temperature in a reaction buffer for initial strand synthesis (5X). First-strand cDNA was generated using random hexamer primers and M-MuLV reverses transcriptase, and RNA was then digested with RNaseH. The DNA polymerase I apparatus was used to create second-strand cDNA. After the 3’ end of the DNA fragment was adenylated, a junction with a hairpin loop structure was added for hybridization. The library fragments were purified using the AMPure XP system (Beckman Coulter, Beverly, USA) to select cDNA fragments with a length of preferably 250–300 BP. After PCR amplification, the PCR products were purified using AMPure XP beads to generate libraries. The library was first measured with a Qubit 2.0 fluorometer, diluted to 1.5 ng/L, and the insert size of the library was determined using an Agilent 2100 Bioanalyzer to ensure that it was with high quality. The effective concentration of the library was accurately determined by RT-PCR once the insert size reached the expected level (the effective concentration of the library was higher than 2 nM) to ensure library quality.

Sequencing and Mapping

Deep sequencing was performed using an Illumina NovaSeq 6000 with paired ends and a 150 bp read length. To ensure accurate data analysis, raw sequences obtained from sequencing were subjected to quality control to obtain good-quality sequences (Clean Reads). The following quality control criteria were applied: removal of reads that contain the junction filtering of low-quality data to ensure data quality, and elimination of reads with over 5% N content (unidentifiable base information). The incoming data was filtered using fastp software based on the aforementioned criteria. Hisat2 v2.2.1 was used to map the Clean Reads to the rat reference genomic RNA from the Ensembl database to the rat reference genome Rattus_norvegicus.mRatBN7.2.109. Finally, featureCounts was used to quantify the bam files to obtain the transcriptome raw counts data.

Data Integration and Analysis

The original counts data from the lab and the counts data of GSE223385 downloaded from the GEO database were combined. The batch effect was removed using the ComBat_seq function of the R package . The combined counts data were transformed into TPM data, and a total of 9 TTS samples and 9 normal samples were used as the training set for analysis. Analysis was conducted using the training set. To confirm the expression of marker genes, the GSE211134 dataset, which includes 20 TTS samples and 9 normal samples, was used. All data were log2-transformed and normalized using the R package. Additionally, the FerrDb database provided the FRGs used in this study. A detailed list of genes is displayed in Table S1.

Differential Expression Analysis

Initially, the expression information for 369 ferroptosis-related genes was collected from the normal and TTS samples in the training set. Next, paired t-tests were used to identify the DE-FRGs in the two distinct samples, and genes with p<0.05 were considered significant.

Functional Enrichment Analysis

GO and KEGG enrichment analyses were conducted on the DE-FRGs using the R package “ClusterProfiler”.

Identification of TTS Diagnostic Genes

The dimensionality of the data was reduced using the Least Absolute Shrinkage and Selection Operator (LASSO) algorithm implemented in the glmnet R package. The LASSO algorithm was used to identify genetic biomarkers of TTS by comparing DE-FRGs in TTS patients and normal samples. Additionally, a Support Vector Machine-Recursive Feature Elimination (SVM-RFE) model was created using the SVM package. The mean false positive rate of 10-fold cross-validation was used to compare the performance of the LASSO and SVM-RFE models. Overlapping genetic markers for TTS were identified using the results of both models. The diagnostic potential of the best genetic biomarker was evaluated by determining the ROC and calculating the area under the curve (AUC), accuracy, sensitivity, and specificity. Furthermore, a logistic regression model based on four marker genes was constructed using the predictive function in the R package glm, to distinguish between different sample types in the training set. The diagnostic performance of the logistic regression model was evaluated using ROC curves.

Gene Set Enrichment Analysis (GSEA)

The R package GSEA (V4.1.2) was utilized for this study. GSEA enrichment analysis was performed separately for the TTS and normal groups to better understand the relevant pathways enriched in TTS and normal samples.

Immune Infiltration Analysis

In this study, the single-sample gene-set enrichment analysis (ssGSEA) algorithm was utilized to estimate the type and relative expression of immune cells that had infiltrated into TTS and normal samples from the training set. The gene set defining the various types of immune-infiltrating cells originates from Charoentong’s comprehensive study, which encompasses a diverse array of human immune cell subtypes. These include activated B cells, activated CD8⁺T cells, activated CD4⁺T cells, and macrophages, among others^16,17 (Table S3).

Quantitative Real-Time PCR

The rat TTS model was reconstructed using the method described in Section 2.1. The GeneJET RNA Purification Kit reagent (Thermo Scientific) was utilized to extract RNA from rat cardiac tissue. The RNA samples were reverse transcribed to cDNA, and the cDNA was by amplified real-time quantitative fluorescence PCR to determine gene expression level. The 2^−ΔΔCT method was employed to analyze the data, and the relative expression of target genes was calculated by the GAPDH ratio. The following is the list of primers: Il1b: Forward: 5’-TAGCAGCTTTCGACAGTGAGG-3’; Reverse: 5’-CTCCACGGGCAAGACATAGG-3’. Enpp2:5’-CAAGGCTGTCTCGGGTCATT-3’; Reverse: 5’-AGAGTCAGACAGCACTGTGG-3’. Etv4:5’-CCGGACAGTGATGAGCAGTT-3’; Reverse: 5’-GGGAGTCATAGGCACTGGAG-3’. GAPDH:5’-GCGGCATTAACTCTGTCCTG-3’; Reverse: 5’-GTTCCCCTTACCTCCAGCTT-3’.

Statistical Analysis

All statistical analyses were carried out using the R software. A Student’s t-test was performed to compare the two groups, while correlation analysis was used to determine the association between DE-FRGs. Statistical significance was set at P <0.05.

Results

Identification of DE-FRGs in the Training Cohort

Out of the 369 FRGs identified in the training set, 37 genes (36 up-regulated and 21 down-regulated) displayed different expression levels between TTS and normal samples. The clustering heatmap demonstrates the expression pattern of these DE-FRGs (Figure 1A). Figure 1B illustrates the interaction between these genes.

Figure 1 (A) DE‐FRGs in normal and TTS samples. (B) The correlation of DE‐FRGs. **P < 0.01, ***P < 0.001.

Functional Analysis of DE-FRGs

Pathway enrichment analysis was conducted by GO and KEGG to identify the signaling pathways associated with DE-FRGs (Table S2). GO analysis discovered that neuronal death and the regulation of neuronal death was closely linked to the DE-FRGs in terms of biological processes, while the outer membrane, organelle outer membrane, and mitochondrial outer membrane were associated with the cellular components. In addition, biological activities linked to the binding of ubiquitin protein ligases, ubiquitin-like protein ligases, and DNA-binding transcriptional repressor activities were found to be related to the DE-FRGs (Figure 2A and B). KEGG pathway enrichment analysis revealed that the DE-FRGs were mostly involved in pathways related to lipid and atherosclerosis, kaposi sarcoma-associated herpesvirus infection, human cytomegalovirus infection, and neurodegeneration multiple disorders (Figure 2C).

Figure 2 (A) GO enrichment analyses. (B) Top 10 enriched GO terms for biological processes (BP), cellular components (CC), and biological functions (BF). (C) Enrichment analysis of the top 30 KEGG pathways.

Four DE-FRGs Were Identified as Diagnostic Genes for TTS

Two different machine learning techniques, LASSO and SVM-RFE, were employed to filter significant DE-FRGs and distinguish TTS from normal samples in the training set. The LASSO regression approach was applied with penalized parameter adjustment by 10-fold cross-validation to select 9 TTS-related characteristics (Figure 3A and B). The SVM-RFE algorithm was then used to screen DE-FRGs and choose the optimal mix of feature genes. Nine genes were selected as the best feature genes with a maximum accuracy of 0.95 and a minimum RMSE of 0.05 (Figure 3C and D). After intersecting the marker genes obtained in the LASSO and SVM-RFE models, Il1b, Enpp2, Pla2g6, and Etv4 were identified as four marker genes for further investigation. Using the R package glm, a logistic regression model was built based on the four marker genes. The ROC curves demonstrated that the model successfully distinguished between normal and TTS samples with an AUC of 1.0 (95% CI: 1.000–1.00). Moreover, ROC curves were plotted for the four marker genes to better understand their ability to differentiate between TTS and normal samples. All genes had an AUC greater than 0.9 (Figure 3E–G). Based on these findings, the logistic regression model is more accurate and specific than individual marker genes in discerning TTS samples from normal samples.

Figure 3 (A and B) By LASSO regression algorithm, with penalty parameter tuning conducted by 10‐fold cross‐validation, was used to select 9 TTS‐related features. (C and D) SVM‐RFE algorithm to filter the 57 DE‐FRGs to identify the optimal combination of feature genes. Finally, 9 genes (maximal accuracy =0.95, minimal RMSE =0.05) were identified as the optimal feature genes. (E) The marker genes obtained from the LASSO and SVM‐RFE models. (F) Logistic regression model to identify the AUC of TTS samples. (G) ROC curves for the 4 marker genes.

Immune Infiltration and GSEA Analysis

The immune infiltration analysis carried out on the training set revealed significant differences in the expression of various immune cells between TTS and normal samples. Mast cells, activated CD4⁺ T cells, eosinophils, myeloid-derived suppressor cells, and activated CD8⁺ T cells were found to be more highly expressed in TTS samples than in normal samples (Figure 4A). The spearman correlation analysis revealed a significant positive correlation between Il1b and CD56bright natural killer cells, CD56dim natural killer cells, mast cells, and MDSC. Conversely, Enpp2 showed a significant negative correlation with activated B cells, CD56bright natural killer cells, CD56dim natural killer cells, macrophages, regulatory T cells, among others (Figure 4B). Furthermore, the results of GSEA analysis indicated that TTS samples were primarily enriched in signaling pathways linked to the cell cycle, cytokine-cytokine receptor interaction, and proteasome, whereas normal samples were mainly enriched in fatty acid metabolism, metabolism of xenobiotics by cytochrome p450, oxidative phosphorylation, and other related signaling pathways (Figure 5A and B).

Figure 4 (A) Exploring differences in immune infiltration between TTS and normal samples using the ssGSEA algorithm. (B) Spearman analysis of the correlation between diagnostic model genes and immune cells. *P < 0.05, **P < 0.01, ***P < 0.001.

Figure 5 (A and B) GSEA analysis was performed on the TTS and normal groups.

Validation of Marker Genes’ Expression

The relevant marker genes were validated in the training set. The expression trends of Il1b, Enpp2, and Etv4 were found to be consistent with the training set. Specifically, Il1b (p<0.001) and Etv4 (p=0.011) expression were observed to be higher in TTS samples than in normal samples, while Enpp2 (p=0.0043) was found to be more highly expressed in normal samples than in TTS samples (Figure 6A–C). Furthermore, we analyzed additional heart tissues from rats in the TTS and normal groups by real-time PCR, and the findings were in line with the previous analysis (Figure 6D–F).

Figure 6 (A–C) The expression of marker genes in the GSE211134 dataset. (D–F) Validating the expression of marker genes using RT-qPCR. **P < 0.01, ***P < 0.001.

Discussion

TTS was initially considered a benign self-limiting disease when it was first reported in 1990. However, recent studies have indicated that it could result in complications such as cardiac arrest, cardiogenic shock, severe mitral regurgitation, ventricular thromboembolism, and life-threatening arrhythmias during the acute phase, with in-hospital mortality rates ranging between 4% and 5%.^18,19 The most popular theory regarding the pathogenesis of TTS suggests that intense mental stimulation or physical stress triggers the activation of the sympathetic nervous system, leading to the secretion of large amounts of catecholamines, which ultimately cause morbidity. Elevated catecholamine levels have been observed in TTS patients both during the acute phase and the recovery period. Moreover, conditions such as pheochromocytoma, subarachnoid hemorrhage, etc. can lead to the development of TTS as well^20,21 Recent evidence also points towards the intricate relationship and interplay between the steroid hormones and catecholamines in the pathogenesis of TTS.²² To establish effective rat model of TTS, we selected the β-agonist isoprenaline as the stimulus as reported before.²³ Initially, we utilized transcriptome sequencing to extract FRGs from the transcriptome data, then the variation in FRGs expression levels between normal and TTS samples was investigated. Two separate machine learning algorithms, LASSO and SVM-RFE, were employed to develop diagnostic models that screen for significant DE-FRGs to differentiate TTS from normal data. In this study, four genes, namely, Enpp2, Pla2g6, Etv4, and Il1b, were identified as differential genes related to ferroptosis. All four genes found to be essential to the diagnostic model have AUC values above 0.9, indicating their high specificity for TTS.Enpp2, also known as autocrine motility factor Autotaxin (ATX), is a hemolytic phospholipase D found in the arterial wall that promotes the production of lysophosphatidic acid (LPAS), leading to atherogenic monocyte adhesion. Enpp2 is involved in a variety of biological functions through the associated G protein-coupled receptor LPAR1-6. Additionally, Enpp2 is a lipid kinase involved in lipid metabolism and has been found to protect cardiac myocytes from elastin-induced ferroptosis. Overexpression of Enpp2 promotes migration and proliferation in H9C2 cells while significantly reducing elastin-induced ferroptosis.²⁴ Pla2g6 is abundantly expressed in human placental trophoblast cells^25,26 and functions as a new ferroptosis regulator in these cells.²⁵ When activated by oxidized phospholipids or other ROS, Pla2g6 can hydrolyze Hp-PE to hemolytic PE and oxidized fatty acids, thus eliminating the ferroptosis signal.²⁷

Furthermore, knocking down Etv4 has been shown to inhibit SLC7A11, a crucial component in ferroptosis regulation, leading to accelerated papillary thyroid cancer growth.²⁸ Inflammatory cytokines such as tumor necrosis factor (TNF), interleukin-1, and interleukin-6 are examples of markers downstream of ferroptosis that can be triggered.²⁹ In our study, we found that the expression of the model gene Il1b was significantly increased in TTS, consistent with the fact that patients with TTS frequently experience an inflammatory response during the acute phase.³⁰ Immune cells play a significant role in maintaining cardiac homeostasis, and recent studies on cardiac injury repair have highlighted the critical role played by immune cells that infiltrate cardiac tissue, such as mast cells, eosinophils, T cells, and B cells, in the preservation of healthy cardiac function.³¹ The iron deficit can prevent T cell proliferation, while iron overload can cause an imbalance between CD4⁺ and CD8⁺ T cell populations. Activated T cells have high iron demand and may limit SLC7A11 expression and impair cystine absorption by tumor cells through IFN-γ production by CD8⁺ T cells, thereby triggering ferroptosis via the JAK1-STAT1 pathway.³² Myeloid-derived suppressor cells (MDSC), which are regulatory immune cells that primarily dampen the immune response,³³ were found to have higher CD4⁺ T cell expression and MDSC infiltration in TTS samples compared to normal samples in our study. These findings suggested that there may be an imbalance in immune cell infiltration in the hearts of rats in the TTS group.

In conclusion, our study effectively developed a four-FRG TTS diagnostic model with potential translational implications for the diagnosis and treatment of TTS, utilizing high-throughput data from our lab and public databases. However, there are certain limitations to this study. Despite reducing batch effects through data integration, confounding variables may still be affecting the results. Additionally, further in vitro and in vivo studies are necessary to confirm the importance of the identified genes.

Conclusion

Our findings support the theoretical importance of ferroptosis in TTS, highlighting Enpp2, Pla2g6, Etv4, and Il1b as potential diagnostic and therapeutic biomarkers for TTS.