Analysis and Validation of Critical Signatures and Immune Cell Infiltration Characteristics in Doxorubicin-Induced Cardiotoxicity by Integrating Bioinformatics and Machine Learning

Abstract

Purpose

Doxorubicin-induced cardiotoxicity (DIC) is a severe side reaction in cancer chemotherapy that greatly impacts the well-being of cancer patients. Currently, there is still an insufficiency of effective and reliable biomarkers in the field of clinical practice for the early detection of DIC. This study aimed to determine and validate the potential diagnostic and predictive values of critical signatures in DIC.

Methods

We obtained high-throughput sequencing data from the GEO database and performed data analysis and visualization using R software, GO, KEGG and Cytoscape. Machine learning methods and weighted gene coexpression network (WGCNA) were used to identify key genes for diagnostic model construction. Receiver operating characteristic (ROC) analysis and a nomogram were used to assess their diagnostic values. A multiregulatory network was built to reveal the possible regulatory relationships of critical signatures. Cell-type identification by estimating relative subsets of RNA transcript (CIBERSORT) analysis was used to investigate differential immune cell infiltration. Additionally, a cell and animal model were constructed to investigate the relationship between the identified genes and DIC.

Results

Among the 3713 differentially expressed genes, three key genes (CSGALNACT1, ZNF296 and FANCB) were identified. A nomogram and ROC curves based on three key genes showed excellent diagnostic predictive performance. The regulatory network analysis showed that the TFs CREB1, EP300, FLI1, FOXA1, MAX, and MAZ modulated three key genes. An analysis of immune cell infiltration indicated that many immune cells (activated NK cells, M0 macrophages, activated dendritic cells and neutrophils) might be related to the progression of DIC. Furthermore, there may be various degrees of correlation between the three critical signatures and immune cells. RT‒qPCR demonstrated that the mRNA expression of CSGALNACT1 and ZNF296 was significantly upregulated, while FANCB was significantly downregulated in DOX-treated cardiomyocytes in vitro and in vivo.

Conclusion

Our study suggested that the differential expression of CSGALNACT1, ZNF296 and FANCB is associated with cardiotoxicity and is also involved in immune cell infiltration in DIC. They might be potential biomarkers for the early occurrence of DIC.

Keywords:

• doxorubicin
• cardiotoxicity
• biomarker
• machine learning
• immune infiltration

Introduction

With the significant advancement in early detection and treatment, the number of cancer survivors has increased,¹ but there have also been multiple adverse outcomes, such as chemotherapy-induced cardiotoxicity.² Doxorubicin (DOX), a member of the anthracycline antibiotic family, is extensively used in the chemotherapy of various cancers, such as hematological malignancies and other tumors.³ Unfortunately, approximately a quarter of cancer survivors experience DOX-induced cardiotoxicity (DIC), which significantly limits its clinical application.⁴ Although DIC has attracted the attention of cardiologists, there are limited methods available for early detection and treatment in clinical practice.

The pathophysiological mechanisms of DIC are complex and not yet fully elucidated, including oxidative stress, mitochondrial dysfunction, inflammatory response, metabolic disorders, and cell death.^5–8 An increasing amount of evidence has indicated that immune cells have crucial functions in cardiac homeostasis and disease.^9,10 Further investigation is required to explore the precise function and molecular mechanism of immune cells in the development of DIC. The rapid developments of multiomics and high-throughput sequencing provide novel tools for identifying early biomarkers and exploring therapeutic targets.¹¹

Human induced pluripotent stem cell-derived cardiomyocytes (hiPSC-CMs) has the ability of unlimited proliferation and differentiation into cardiomyocytes, so it is widely used in the study and treatment of cardiac disease.^12,13 Recently, Burridge et al indicated that the viability of using hiPSC-CMs as a platform for determining and validating the genetic foundation and molecular mechanisms underlying DIC.¹⁴ In addition, many studies have also used hiPSC-CMs as the cell source of DIC.^15–17 On these bases, we downloaded high-throughput gene expression data from NCBI GEO database.

In the current work, we explored key genes of DIC using comprehensive bioinformatics approaches and machine-learning strategies. Furthermore, we constructed a multiregulatory network to uncover the possible regulatory connections of key genes. Similarly, machine learning was carried out to construct a diagnostic nomogram model for DIC prediction on the basis of the three key genes (ZNF296, FANCB and CSGALNACT1) that were discovered in DIC pathogenic genes. Therefore, the study of its related molecular mechanism may be helpful for its further clinical application in DIC. We attempted to explore the connection between key genes and immunity. Finally, our study validated the differential expression of the key genes by constructing the model treated with doxorubicin in vitro and in vivo. This study collectively identified a stable and achievable biomarker in DIC, which could potentially serve as a target for diagnosing and preventing cancer survivors who have undergone DOX treatment.

Materials and Methods

Data Acquisition and Homogenization

Our information was gathered by accessing the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/), which is a public repository of high-throughput gene expression data, chips, and microarrays.¹⁸ High-throughput gene expression data were downloaded from the NCBI GEO database, including GSE76314, GSE157282 and GSE181517. GPL11154 served as the complementary high-throughput detection platform for GSE76314; GSE157282 was based on platform GPL24676. GSE181517 was based on platform GPL18573. The GSE76314 dataset,¹⁴ which describes the differences in cell-derived cardiomyocytes (hiPSC-CMs), was used as the training set (3 doxorubicin samples vs 3 control samples). The GSE157282 dataset¹⁹ [doxorubicin samples vs 3 controls samples], which describes the differences in the hiPSC-CMs, was also used as the training set. GSE181517²⁰ was used as the validation set to explore the dynamic changes in the critical signature, which included 3 doxorubicin samples and 3 control samples. The scale function in R version 4.2.0 software was used to combine GSE76314 and GSE157282 into one dataset and calibrate the batch using the SVA package.²¹ For uniform manifold approximation and projection (UMAP), after the removal of the batch effect, the samples among various datasets are clustered and intertwined with each other, indicating that the batch effect is better removed. Principal component analysis (PCA) was used to verify the reproducibility of the data, and the R package ggord was used to construct the PCA plots.²²

Analysis of Differentially Expressed Genes

In this study, we performed differential expression analysis using the R Bioconductor package “Limma” (V3.46.0) to identify differentially expressed genes.²³ Differentially expressed genes (DEGs) were defined as those with adjusted P values and genes with adjusted P values < 0.05 and |Log2-fold-change (log2FC)|≥ 1 after adjusting the P values using Benjamini‒Hochberg’s false discovery rate (FDR). The R software ggplot2 V3.3.5 package was used to create volcano plots, and the R software Heatmap V1.0.12 package was utilized to generate heatmaps for the top 100 DEGs from each dataset.

Functional Enrichment Analysis of DEGs

A functional enrichment analysis of DEGs was performed to clarify their role using the R software clusterProfiler and Goplot package. These analyses were discerningly marked by significance, indicated by the adjusted P value < 0.05. A comprehensive enrichment analysis of Kyoto Encyclopedia of Genes and Genomes (KEGG) and gene ontology (GO) was conducted and visualized, including biological process (BP), cellular component (CC), and molecular function (MF) for all DEGs.

Screening and Validation of Critical Gene Signatures

In this study, random forests (RF),²⁴ least absolute shrinkage and selection operator (LASSO) logistic regression,²⁵ and weighted gene coexpression network analysis (WGCNA)²⁶ were utilized to screen novel and pivotal signatures associated with DIC. The RF model, LASSO analysis, and WGCNA were executed using the R software random Forest package, R software glmnet package and R software WGCNA package, respectively. Following the individual algorithmic analyses, the points of intersection among the three classification models were identified, setting the stage for the next step of the analysis. To assess the diagnostic efficacy of these three diagnostic signatures in doxorubicin-induced cardiotoxicity, the area under the curve (AUC) and the ROC curves were established using R software’s pROC package. In adherence to statistical standards, a threshold of two-sided P < 0.05 was deemed to be statistically significant. To verify the expression levels of key markers, their differences between the doxorubicin and control groups were assessed in the training and validation sets.

Multifactor Regulatory Network Analysis of Key Genes

Guardians of gene expression include long noncoding RNAs (lncRNAs), microRNAs (miRNAs), and transcription factors (TFs) in transcription and posttranscription prototypes. To demonstrate the regulatory networks of key signatures, lncRNAs, mRNAs, TFs, and miRNAs were predicted by RNAInter,²⁷ hTFtarget,²⁸ and mirDIP.²⁹ Then, we used Cytoscape software³⁰ to visualize the multifactor regulatory network.

Analysis of Immune Cell Infiltration on Critical Gene Signatures

Prior to CIBERSORT analysis, principal component analysis (PCA) of the gene expression profile was plotted and visualized by the R software factoextra package. Analysis of 22 different types of immune cells using CIBERSORT demonstrated the differential expression of cells between the doxorubicin group and the control group using a box diagram (P value < 0.05 was deemed to be statistically significant). Afterwards, Spearman correlation analysis between immune cells and immune cells, as well as between key signatures and immune cells, was demonstrated by a heatmap through the R software corrplot V0.92 package.

Cell Line Culture and Treatment

AC16 immortalized human cardiomyocyte line was procured from MeisenCTCC (Zhejiang, China) and cultured in Dulbecco’s modified Eagle’s medium (DMEM) containing 10% fetal bovine serum (FBS) 1% and penicillin/streptomycin and grown in a CO₂ incubator maintained at atmospheric oxygen levels and 5% CO₂. When the AC16 cell density reached 80–90%, the cells were exposed to 1 µM doxorubicin (Sigma‒Aldrich) or vehicle (PBS). Whole cells were collected after 24 h of treatment.

Animal Experiments

The Establishment of Animal Model

Twenty male wild-type C57BL/6 mice aged 8 weeks (18–22g) were purchased from the SiPeiFu (Beijing) Biotechnology Co., Ltd., Beijing City, and they were randomly divided into two groups: normal control (NC) group and DOX treatment (DOX) group. In NC group, mice were i.p.injected with the equivalent volume of saline two times (on days 1 and 8, respectively); in DOX group, mice were intraperitoneally injected with DOX on day 1 (10 mg/kg) and day 8 (10 mg/kg), and acute DIC model was established in vivo (cumulative dose of 20 mg/kg).³¹ At the end of the experiment, each group of mice was anesthetized and killed, and their hearts were taken for subsequent experiments. The experimental procedures were implemented in accordance with the China Animal Welfare Legislation, reviewed and approved by the Qingdao University Committee on Ethics in the Care and Use of Laboratory.

Cardiac Functions Examination

At the end of the treatment, the mice heart function was examined through echocardiography Imaging System (VINNO 6 LAB, VINNO). Cardiac ventricular structure and systolic function were examined by two-dimensional imaging motion mode to record parasternal long-axis fields. Parameters included the left ventricular (LV) anterior wall thickness in the end-diastolic/end- systolic stages, the LV posterior wall thickness, the LV internal diameter; LV ejection fraction, and LV fractional shortening.

The Hematoxylin & Eosin Staining of Heart Tissues

Hearts of mice were fixed in 10% formalin, and then treated with dehydration, embedding and sectioning. Hematoxylin & eosin staining was performed to examine cardiac histological morphology.

RNA Extraction and Quantitative Real-Time PCR for Cell and Animal Verification Experiments

The TRIzol method was used to extract total RNA from AC16 cells and heart tissue samples, and spectrophotometry was used to measure the concentration and purity. According to the manufacturer’s instructions, equivalent amounts (1 μg) of purified RNA were used as a template to synthesize cDNA using ABScrlpt III RT Master Mix for qPCR with a gDNA Remover Kit (Cat# RK20429, ABclonal, Wuhan, China). The manufacturer’s instructions were followed to perform quantitative PCR using a QuantStudio 3 (Bio-Rad) and 2X Universal SYBR Green Fast qPCR Mix (Cat#RK21203, ABclonal, Wuhan, China). The relative expression levels were determined by the 2^–ΔΔCT method and are presented relative to GAPDH mRNA. Primer sequences are provided in Table S1.

Statistical Analysis

GraphPad Prism 9.0 was utilized for conducting the statistical analysis. The unpaired Student’s t-test was used to compare the DOX group (DOX) and control group (NC). The assays were conducted at least three times or more. A two-tailed p value <0.05 was considered statistically significant.

Results

Identification of DEGs

Figure 1 clarifies the depicted flow chart of the present work. Both the GSE76314 and GSE157282 datasets were merged (Figure 2A) for the following analysis. The utilization of UMAP was performed to visually represent the clustering of read counts and identify batch effects (Figure 2B). Next, we merged these two datasets into one dataset and rectified the batch using the SVA software package. Figure 2C shows that the interbatch variations are effectively removed after data normalization. Finally, a normalized gene expression matrix file containing 12 samples (6 doxorubicin samples and 6 control samples) was obtained. These DEGs are distributed in the volcano map (Figure 2D). The first 100 temperatures are shown in the heatmap (Figure 2E). In the results of differential expression analysis, a total of 3713 genes were identified as DEGs, including 1925 upregulated genes and 1788 downregulated genes (Figure 2F).

Figure 1 Flow chart of this study.

Figure 2 Identification of differentially expressed genes (DEGs). (A) The intersection of genes between GSE76314 and GSE157282. (B) UMAP before the batch correction of two datasets. (C) UMAP after the batch correction of two datasets. (D) Volcano plots of DEGs. Nodes in red represent upregulated genes, nodes in green represent downregulated genes, and black dots represent no significantly changed genes. (E) Heatmaps of DEGs. The legend on the top right indicates the log fold change of the genes. The horizontal axis represents each sample, and the vertical axis represents each gene. Blue and red colors represent low and high expression values, respectively. (F) Number of DEGs.

GO and KEGG Pathway Analysis of DEGs

Subsequently, we performed enrichment analyses for overlapping DEGs in KEGG and GO. Among the KEGG pathways, the pathways in cancer, erpes simplex virus 1 infection, rap1 signaling pathway, alcoholism, and cGMP-PKG signaling pathway were considered the pathways with the highest degree of enriched pathways (Figure 3A). In significantly enriched BP, these included the positive regulation of nucleobase containing compound metabolic process, positive regulation of biosynthetic process, cell population proliferation, cell cycle, and apoptotic process (Figure 3B). Subsequently, chromosome, chromatin, microtubule cytoskeleton, transcription regulator complex, and chromosomal region were the most enriched terms in CC (Figure 3C). For MF, the top 5 enriched items were transcription regulator activity, sequence-specific DNA binding, ribonucleotide binding, DNA binding transcription factor activity, and adenyl nucleotide binding (Figure 3D).

Figure 3 Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of differentially expressed genes (DEGs). (A) Bubble chart shows KEGG-enriched items of DEGs. (B–D) Bubble charts show GO-enriched items of DEGs in three functional groups: biological processes (BP, B), cell composition (CC, C), and molecular function (MF, D). The x-axis labels represent gene ratios, and the y-axis labels represent GO terms. The size of the circle represents the gene count. Different colors of circles represent different adjusted P values.

Screening of Critical Signatures

Here, three key genes were screened using the RF, LASSO, and WGCNA machine learning algorithms. Using the WGCNA algorithm (Figure S1), we confirmed 11 extraordinary coexpression modules (Figure 4A). Figure 4B clarifies that there is a correlation between multiple modules in DIC. Given that the turquoise and dark red modules are most significantly associated with DIC, we screened genes in the turquoise and dark red modules for follow-up exploration and successfully identified 3370 genes (Figure 4C–E, Table S2). Moreover, the results showed that 15 genes were identified by the LASSO algorithm (Figure 4F and G, Table S2), and 90 genes were determined with the RF algorithm (Figure 4H–J, Table S2). Subsequently, a Venn diagram shows that FANCB, ZNF296 and CSGALNACT1 obtained by the three algorithms are overlapping genes (Figure 4K).

Figure 4 Screening of critical signatures via multiple machine-learning algorithms. (A) Clustering dendrogram of differentially expressed genes related to DIC, with dissimilarity based on topological overlap, together with assigned module colors. (B) Module feature vector clustering. (C) Module–trait associations. The gene significance for doxorubicin cardiotoxicity in the turquoise and dark red modules, and one dot represents one gene in the turquoise (D) and dark red modules (E). (F) LASSO coefficient profile of the 12 genes, and different colors represent different genes. (G) Selection of the optimal parameter (lambda) in the LASSO model and generation of a coefficient profile plot. (H) Variable importance, as measured by the mean decrease in accuracy (left panel) or the Gini coefficient (right panel), is computed using the OOB error. Genes are shown in descending order of importance. Distribution of the out-of-band (OOB) error rate at various values of mtry (I) and trees (J). (K) Venn diagram showing the intersection of critical signatures obtained by the three strategies.

Verification of FANCB, ZNF296 and CSGALNACT1 in DIC

To evaluate the possible significance of the key genes in DIC, we created ROC curves. First, we verified the expression of FANCB, ZNF296 and CSGALNACT1 in DIC. The results showed that compared with the control group, the expression of ZNF296 and CSGALNACT1 in the training and validation sets was significantly upregulated, and the expression of FANCB was significantly downregulated (all P < 0.01, Figure 5A and B). Second, the AUCs of the three key genes in the training set were 1 (Figure 5C), indicating that they were accurate and reliable in terms of predictive value. In the validation set of GSE181517, the AUCs for FANCB, ZNF296 and CSGALNACT1 were also all 1 (Figure 5D). Finally, we created nomogram models to diagnose DIC using the hallmark genes FANCB, ZNF296 and CSGALNACT1 (Figure 5E).

Figure 5 The diagnostic power of FANCB, ZNF296 and CSGALNACT1 in DIC by ROC curve. (A) The expression of FANCB, ZNF296 and CSGALNACT1 in the training set. (B) The expression of FANCB, ZNF296 and CSGALNACT1 in the validation set. (C) The ROC curve of FANCB, ZNF296 and CSGALNACT1 in the training set. (D) The ROC curve of FANCB, ZNF296 and CSGALNACT1 in the validation set. (E) A nomogram was used to predict the occurrence of DIC. (*P<0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001, - P ≥ 0.05).

Construction of a Multifactor Regulatory Network Based on Key Genes

Next, we explored the potential regulatory mechanisms by studying the regulatory networks involving lncRNAs, miRNAs and TFs and conducted a comprehensive visual analysis (Figure 6). Here, to construct a multifactor regulatory network, interaction pairs involving miRNAs, lncRNAs and TFs with three key genes were extracted. In the multifactor regulatory network, FANCB was managed by 5 miRNAs and 1 lncRNA, including hsa-miR-4695-5p, hsa-miR-3612, hsa-miR-2113, hsa-miR-7109-3p, hsa-miR-135a-2-3p, and CISTR, and the possible TFs for FANCB were ATF3, CEBPA, CEBPB, CREB1, E2F1, EP300, ETS1, FLI1, FOXA1, GABPA, HDAC1, HDAC2, KLF5, MAX, and MAZ. Additionally, the ZNF296 network included 2 lncRNAs, 7 miRNAs, and 66 TFs, and the CSGALNACT1 network contained 167 miRNAs and 25 TFs. Obviously, the results of our work suggest that the TFs CREB1, EP300, FLI1, FOXA1, MAX, and MAZ modulate all 3 marker genes.

Figure 6 The multifactor regulatory network based on FANCB, ZNF296 and CSGALNACT1. FANCB, FA complementation group B; ZNF296, zinc finger protein 296; CSGALNACT1, chondroitin sulfate N-acetylgalactosaminyltransferase 1.

Analysis of Immune Cell Infiltration and Correlation Analysis Between Key Genes

Recent evidence has shown that immune cells play a critical role in the pathophysiological process of DIC.^6,32,33 Therefore, we investigated the link between key features in DIC and immune infiltration. The swatches of different groups are well divided in the training set in PCA (Figure 7A). To study 22 types of immune cell phenotypes, we utilized the CIBERSORT algorithm to analyze the training set. Compared with the control samples, DIC samples exhibited a greater percentage of activated NK cells, M0 macrophages, activated dendritic cells and neutrophils (all P < 0.05). However, the proportions of memory B cells, CD8 T cells, resting memory CD4 T cells, and resting NK cells in DIC were relatively lower than those in the control (all P < 0.05) (Figure 7B). Twenty immune cell types were included for further investigation after two immune cell types with undetectable levels were excluded. It can be seen from the related heatmaps (Figure 7C) that naive B cells and CD8 T cells, naive B cells and resting memory CD4 T cells, CD8 T cells and neutrophils, resting NK cells and activated M0 macrophages displayed the most significant negative correlations (P < 0.0001). Naive B cells and neutrophils, activated mast cells and monocytes, and M0 macrophages and neutrophils exhibited the most significant positive correlations (P < 0.0001). Our aim was to detect the link between key genes and 20 types of immune cell phenotypes in doxorubicin cardiotoxicity. According to the correlation analysis results, FANCB displayed the most positive correlations with resting NK cells and CD8 T cells [r = 0.70, -log10 (p value) = 1.94] while demonstrating negative correlations with activated dendritic cells [r = −0.76, -log10 (p value) = 2.35]. ZNF296 displayed the most positive correlations with activated dendritic cells [r = 0.76, -log10 (p value) = 2.40] and revealed negative correlations with resting NK cells [r = −0.80, -log10 (p value) =2.72]. CSGALNACT1 displayed the most positive correlations with M0 macrophages [r = 0.92, -log10(p value) = 4.60] and revealed negative correlations with resting NK cells [r = −0.79, -log10(p value) =2.65] (Figure 7D).

Figure 7 Immune cell infiltration analysis and relationships between critical signatures and immune cells in DIC. (A) Principal component analysis (PCA) cluster plot of gene expression profiles between doxorubicin samples and control samples in the training set. (B) Box plot of the proportion of 22 types of immune cells. (C) Heatmap of correlations in 20 types of immune cells. The colored squares represent the strength of the correlation; red represents a positive correlation, and blue represents a negative correlation. Darker color implies stronger association. (D) Correlations between FANCB, ZNF296, CSGALNACT1, and infiltrating immune cells. (*P<0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001, - P ≥ 0.05).

Verification of FANCB, ZNF296 and CSGALNACT1 in AC16 Human Cardiomyocyte-Like Cells

To verify whether the mRNA expression of FANCB, ZNF296 and CSGALNACT1 was consistent with the results of the database analysis, we constructed a doxorubicin-induced injury model in vitro.³⁴ The results showed that FANCB was expressed at a low level (Figure 8A), while ZNF296 and CSGALNACT1 were expressed at a high level in the DOX group (Figure 8B and C).

Figure 8 The mRNA expression of FANCB (A), ZNF296 (B) and CSGALNACT1 (C) in AC16 human cardiomyocyte cells by RT‒qPCR. (**P < 0.01, ***P < 0.001).

The Verification Experiments in vivo

Firstly, mouse model of DIC is illustrated in Figure 9A, and the heart weight to tibia length ratio (HW/TL) was significantly reduced in the DOX group (Figure 9B). The echocardiography results demonstrate a downward trend of cardiac systolic and diastolic function in the DOX group (Figure 9C), as evidenced by significantly decreased ejection fraction (Figure 9D), decreased fractional shortening (Figure 9E). HE staining illustrated that DOX treatment exacerbated cardiomyocyte tissue disorders compared with control group (Figure 9F). These results indicate that we successfully established a DIC model. Finally, the mRNA levels of FANCB was decreased in the hearts of DOX treated group relative to the control group (Figure 9G), while ZNF296 and CSGALNACT1 is contrary to the FANCB (Figure 9H and I).

Figure 9 The mRNA expression of FANCB, ZNF296, and CSGALNACT1 in DOX-induced cardiotoxicity in vivo. (A) Male C57BL/6 mice were i.p.injected with an accumulation does of 20 mg/kg DOX or the equivalent volume of saline two times (on Days 1 and 8, respectively). (B) the HW/TL ratio was measured. (C) The mice were tested for cardiac function by echocardiography. (D and E) The LV trace EF% and LV trace FS% were calculated. (F) HE staining of heart tissue in Control group and DOX group. (G) FANCB mRNA expression was detected by RT-qPCR, (H–I) ZNF296, and CSGALNACT1 mRNA expression was detected by RT-qPCR. (*P<0.05, **P < 0.01, ****P < 0.0001).

Discussion

Doxorubicin prolongs the survival time of tumor patients, but it also leads to severe side effects, especially cardiotoxicity, which significantly affects the effectiveness of anticancer treatment.^35,36 Hence, it is crucial to recognize secure and effective biomarkers to detect DIC. Currently, with the rapid development of high-throughput sequencing and multiomics technologies,^37,38 some new methods have been used to detect and diagnose DIC.¹¹

In this study, 3713 DEGs were discovered using a comprehensive bioinformatics analysis between DIC and the controls. The results indicated that the significantly enriched terms included positive regulation of nucleobase-containing compound metabolic process in BP, chromosome in CC, and transcription regulator activity in MF. Of note, multiple algorithms were applied to identify key genes in DIC, including WGCNA, RF and LASSO. Ultimately, CSGALNACT1, ZNF296 and FANCB were determined to be key genes. Additionally, through ROC curves and nomogram analysis, it was shown that the expression levels of CSGALNACT1, ZNF296 and FANCB could accurately distinguish DIC from controls. Similarly, we validated the stability of CSGALNACT1, ZNF296 and FANCB expression levels through distinct validation datasets. From the results, CSGALNACT1, ZNF296 and FANCB may be potential factors for DIC, which suggests that they have vital value for clinical application in DIC. Additionally, our regulatory network showed that the TFs CREB1, EP300, FLI1, FOXA1, MAX, and MAZ modulated all 3 marker genes. However, the transcriptional regulation of CSGALNACT1, ZNF296 and FANCB by CREB1, EP300, FLI1, FOXA1, MAX and MAZ needs further study.

Zinc finger protein 296 (ZNF296 or Zfp296) is a reprogramming factor composed of 445 amino acids that belongs to the C2H2 zinc finger protein family and is involved in cell proliferation, survival, differentiation, and carcinogenesis.³⁹ ZNF296 is expressed in mouse testes, bone marrow, and embryos.⁴⁰ Previous evidence suggests that ZNF296 is an effective biomarker for breast cancer detection.⁴¹ Furthermore, recent studies have found that ZNF296 is a regulator of H3K9me3 and is closely associated with the pluripotency of embryonic stem cells.⁴² Another study indicated that ZNF296 was associated with the formation of induced pluripotent stem (iPS) cells.⁴³ Interestingly, a microarray analysis based on transcriptomic signatures showed that ZNF296 was upregulated in tuberculosis and could be used as a diagnostic biomarker.⁴⁴ In the present study, we also found that ZNF296 was significantly upregulated in the DOX-treated group, which was validated through platforms and across cell experiments. Although the role of ZNF296 in DIC has not been clearly discussed, it may be a potential diagnostic biomarker.

Furthermore, CSGALNCT-1, also named chondroitin sulfate N-acetylgalactosaminyltransferase 1, plays a vital role as a glycosyltransferase in the biosynthesis of chondroitin sulfate. It is also involved in important biological processes, such as cartilage development and polysaccharide metabolism.⁴⁵ Previous microarray studies have indicated that CSGALNACT1 is a prognostic biomarker in multiple myeloma.⁴⁶ According to a recent study, CSGALNACT1 was found to have a notable upregulation in lung-derived mesenchymal stem cells (LMSCs) from emphysema patients, which can regulate extracellular matrix (ECM), but the specific function remains unclear.⁴⁷ Significantly, chondroitin sulfate (CS) is a natural macromolecule catalyzed by CSGALNACT1 and a major component of the extracellular matrix.⁴⁸ Another study published recently showed that abnormal accumulation of chondroitin sulfate glycosaminoglycans led to cardiac fibrosis.⁴⁹ Interestingly, our study showed that CSGALNACT1 was expressively upregulated in cardiomyocytes in the doxorubicin-treated group. Therefore, we speculate that CSGALNACT1 may be involved in DIC, although the specific mechanisms require further investigation. In conclusion, CSGALNACT1 is considered a suitable biomarker for DIC.

FANCB, fully defined as Fanconi anemia complementation group B and known as FAAP95, is an essential constituent of the Fanconi anemia (FA) core complex. It has a vital function in hematopoiesis and the development of germ cells.^50,51 Some studies have indicated that the FA pathway serves as an indispensable mechanism for the DNA damage response.⁵² Further study has shown that FANCB is known as a DNA repair gene.⁵³ Previous studies have concluded that DNA damage is one of the molecular mechanisms in DIC. Interestingly, our study found that the expression of FANCB in cardiomyocytes treated with doxorubicin was significantly reduced, which was consistent with the verification of the cell experiment. Therefore, we speculate that FANCB is directly or indirectly involved in DIC via DNA damage/repair. Although the exact mechanism remains unclear, this study provides a novel perspective for understanding doxorubicin-induced cardiomyocyte injury.

Noncoding RNAs (ncRNAs) widely participate in cell physiological and pathological processes by regulating the epigenetic modification of gene expression, which has attracted great attention.^54,55 Increasing evidence has reported the significant role of ncRNA in the DIC.⁵⁶ Therefore, we further predicted upstream ncRNAs targeting three key genes, ZNF296, FANCB, and CSGALNACT1, and ultimately constructed a gene regulatory network including microRNAs and lncRNAs. Additionally, we also investigated the regulatory network between transcription factors and three core genes. Notably, the transcription factors CREB1, EP300, FLI1, FOXA1, MAX, and MAZ were simultaneously regulated by three key genes.

In recent years, a growing body of research has highlighted the potential impact of immune cells on DIC, including neutrophils,^57,58 macrophages,^59,60 and NK cells.⁶¹ In our study, significant differences were observed between the DOX-treated group and the control group in terms of immune cell infiltration. Specifically, the peaks of M0 macrophages, activated NK cells, activated dendritic cells, and neutrophils were higher, while the proportions of memory B cells, CD8 T cells, resting CD4 T cells, and resting NK cells were lower. Previous studies have shown that 24 hours after doxorubicin treatment, neutrophils significantly infiltrate the heart, accompanied by a decrease in cardiac function, vascular structure damage, and increased deposition of vascular collagen, leading to fibrosis,⁵⁸ which is consistent with our study findings. Notably, macrophages, a type of innate immune system cell, were also found to be upregulated in our study. Some studies found that in an acute doxorubicin-induced cardiotoxicity mouse model, doxorubicin increased the proportion of M1 macrophages while suppressing M2 macrophages.⁶² A macrophage profile tracking study based on CX3CR1 found that M1 macrophages may be the major group during the early stage of cardiac injury.⁶³ Although studies have suggested that immune cells play an important role in DIC, the specific regulatory mechanisms remain unclear.

Therefore, we attempted to reveal the potential function of immune cells by establishing a connection between key genes and immune cells. In our research, CSGALNACT1, ZNF296 and FANCB showed different degrees of correlation with immune cells, including resting NK cells and activated dendritic cells. In a prior study, it was demonstrated that CSGALNACT1 was linked to the activation of lymphocytes and the production of inflammatory cytokines (IL-6 and IFN-r cytokine) in experimental autoimmune encephalomyelitis (EAE).⁶⁴ In our data, CSGALNACT1 displayed the most positive correlations with M0 macrophages and showed negative correlations with resting NK cells. Our findings provide novel insights into the immune mechanisms of CSGALNACT1, ZNF296, and FANCB during DIC.

There are some limitations in our study. First, although we have integrated multiple datasets, the sample size was still relatively small. Second, due to the lack of relevant clinical data, we did not perform prognostic analysis. Third, we validated the mRNA expression of critical key genes in DOX-treated cardiomyocytes, but further validation, such as Western blotting and immunofluorescence staining, is necessary. To improve these deficiencies, we will conduct further clinical studies to confirm the findings with more detailed clinical information and a wide range of samples.

Conclusion

In summary, we concluded that (1) three key genes, CSGALNACT1, ZNF296, and FANCB, were identified that might be probable diagnostic biomarkers in DIC and may contribute to immunity regulation simultaneously; (2) immune cell infiltration was closely related to DIC; and (3) the inflammatory response signaling pathway, biosynthetic process, and apoptotic process may be the key mechanisms in DIC. This study provides a new strategy for screening and evaluating biomarkers in clinical practice while also laying the foundation for the underlying role of immune cell responses in DIC.

Abbreviations

DIC, Doxorubicin induced cardiotoxicity; mRNA, Message RNA; TCGA, The Cancer Genome Atlas; GEO, Gene Expression Omnibus; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; AUC, Area under the curve; DEGs, differentially expressed genes; lncRNA, Long noncoding RNA; miRNA, microRNA; TF, Transcription factor; ROC, Receiver operating characteristic; RF, Random forests; LASSO, Least absolute shrinkage and selection operator; WGCNA, Weighted gene coexpression network analysis; FANCB, FA complementation group B; ZNF296, Zinc finger protein 296; CSGALNACT1, Chondroitin sulfate N-acetylgalactosaminyltransferase 1; DMEM, Dulbecco’s modified Eagle’s medium.

Data Sharing Statement

The following accession numbers were included in this study that were downloaded and analyzed in the public datasets and could be accessed in the GEO data repository (https://www.ncbi.nlm.nih.gov/geo/): GSE76314, GSE157282, and GSE181517. Further enquiries can be directed to the corresponding author.

Ethics Committee Approval

The animal experiments were approved by the Ethics Committee of Qingdao University.

Author Contributions

All authors made a significant contribution to the work reported, whether that is in the conception, study design, execution, acquisition of data, analysis and interpretation, or in all these areas; took part in drafting, revising or critically reviewing the article; gave final approval of the version to be published; have agreed on the journal to which the article has been submitted; and agree to be accountable for all aspects of the work.

Disclosure

All the authors declared that they had no competing interests for this work.

Acknowledgments

Thank you for the help of the members from the Department of Cardiology at the Affiliated Hospital of Qingdao University.

Additional information

Funding

This work was supported by The National Natural Science Foundation of China (no. 82172574, no. 81871231), The Natural Science Foundation of Shandong Province (no. ZR2020MH016, no. ZR202209280042), Research Planning Project of Shandong Higher Medical Education Research Center (no. YJKT202171), and Shandong Taishan Scholars Young Experts Program (No. tsqn202103056).