This article is mentioned in:

Introduction

Heart failure (HF) is a long-standing public health issue, causing high morbidity and mortality worldwide (1). More than 30 million individuals are currently suffering from HF, and this number is expected to continue rising owing to the global aging population (2). HF is characterized by the inability of the heart to pump sufficient blood for the metabolic demands of the body, resulting in dyspnea, fatigue, poor exercise tolerance and fluid retention (3,4). Therefore, HF poses a considerable economic burden to global healthcare systems (1,5-7). The development of HF is usually accompanied by the synchronous reprogramming of gene expression, associated with modifications in multiple genes and signaling pathways involved in HF pathogenesis (8-11). Dilated cardiomyopathy (DCM), one of the leading factors associated with HF, is characterized by left ventricular dilation combined with systolic dysfunction (12-15). However, the overall picture of myocardial gene co-expression signatures in DCM-induced HF remains unclear. Therefore, detecting key genes and pathways in the pathogenesis of DCM-induced HF is necessary and meaningful, and may perform an important role in the prevention, diagnosis and treatment of DCM-induced HF.

With the development of second-generation sequencing technologies, several gene chip techniques have been used in both research and clinical settings for the treatment of cardiovascular diseases (16-19). Clarke et al (16) found that two lipoprotein variants were associated with an increased risk of coronary disease using a novel gene chip containing single-nucleotide polymorphisms. Kuehl et al (17) identified distinctive microRNAs to assess the risk of virus persistence and progressive clinical deterioration in the course of enterovirus cardiomyopathy using microRNAs gene chips. Weighted gene co-expression network analysis (WGCNA) has been introduced as a novel and powerful systems biology method designed for constructing a co-expression network among identified genes (20). WGCNA divides genes into separate modules according to gene expression patterns that connect with the complex changes of the clinical phenotype. Compared with other bioinformatics analysis methods, this technique offers a more convenient and effective way to obtain key modules that are more closely related to clinical traits; subsequently, hub genes and key pathways can be identified within these modules (20,21). WGCNA has been widely used to identify novel biomarkers and pathways in a number of cardiovascular diseases, including human atrial fibrillation, ischemic cardiomyopathy and acute myocardial infarction (22-24).

In the present study, WGCNA was performed to identify key genes and biological pathways associated with DCM-induced HF. The potential function of the key module was explored using Gene Ontology (GO) annotation, Gene Set Enrichment Analysis (GSEA) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses with the clusterProfiler package in R. The diagnostic efficacy of key genes was evaluated using a receiver operating characteristic (ROC) curve. The key genes were also verified in the GSE116250 dataset and in vitro experiments. Furthermore, the analysis of infiltrated immune cells in heart tissues with DCM-induced HF was performed using the CIBERSORT method.

Materials and methods

Data collection and processing

A flow chart of the present study is presented in Fig. 1. Microarray datasets GSE79962 and GSE116250 were extracted from the GEO database (https://www.ncbi.nlm.nih.gov/geo). The GSE79962 dataset, comprised of 9 patients with DCM-induced HF and 11 non-failing donors, was used to explore novel biomarkers (25). The GSE116250 dataset comprised 37 patients with DCM-induced HF and 14 control patients and was served as the independent testing set to validate the expressions and diagnostic efficacy of key genes (26,27). The selection criteria for these datasets were: i) Microarray dataset had to be available for both patients with DCM-induced HF and healthy controls; and ii) an annotation file for the platform was required for each dataset. All data analyses and processing were conducted using R software.

Figure 1 Flow chart of the study. In this study, two DCM-induced HF datasets were collected. WGCNA and differentially expressed genes analysis were used to identify key genes in the GSE79962 dataset. The relative fractions of infiltrated immune cells were evaluated the CIBERSORT algorithm. The relationship between key genes and immune cells was calculated by Spearman's correlation analysis. The potential function of the key module was explored using GO, GSEA and KEGG analysis. Moreover, key genes were validated in the GSE116250 dataset and in vitro experiments. The diagnostic efficacy of key genes was evaluated using a ROC curve in the GSE116250 dataset. DCM, dilated cardiomyopathy; GO, gene ontology; GSEA, Gene Set Enrichment Analysis; HF, heart failure; KEGG, Kyoto encyclopedia of genes and genomes; ROC, receiver operating characteristic; WGCNA, weighted gene co-expression network analysis.

Co-expression network construction

Co-expression networks were established for 4,343 genes (the top 25% of rank genes with the largest variance) in the GSE79962 dataset using the WGCNA package v1.70-3 (https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA). The soft threshold (β) was set to 9 to construct a co-expression network that conformed to the scale-free distribution when the degree of independence was 0.9. Next, an adjacency matrix was constructed by raising the correlation matrix to the power of 9, and then a topological overlap matrix (TOM) was used to measure similarity based on the adjacency matrix. Genes were hierarchically clustered and visualized in a dendrogram according to the dissimilarity TOM. Hierarchical clustering was performed to obtain modules using a dynamic tree cutting algorithm with a minimum module size of 30 for the gene dendrogram. The modules with similar expression profile were merged, and 24 modules were subsequently obtained.

Functional enrichment of the key modules

Gene significance (GS) is defined as the absolute value of the correlation between between the gene expression and clinical traits. Module significance (MS) is the average GS in a specific module, which represents the correlation between the module and clinical traits (20). The module with the highest MS value was considered the key module most relevant to DCM-induced HF. To explore the potential mechanisms and functions of the target modules, GO functional term and KEGG pathway enrichment analyses were performed using the R package ‘clusterprofiler’ (v 4.0.5) (28).

GSEA

The related biological pathways in DCM-induced HF from the GSE79962 dataset were further analyzed using GSEA in the R package ‘clusterprofiler’ (v 4.0.5).

Differentially expressed genes (DEGs)

The ‘limma’ R package was used to identify DEGs between healthy individuals and patients with DCM-induced HF in the GSE79962 dataset. An adjusted P<0.05 and log(fold change) >1.5 were set as the cut-off values for DEG screening.

Identification of key genes

Key genes were mined according to the threshold value of module membership (MM), GS and DEGs analysis. MM refers to the Pearson's correlation coefficient between genes and the module eigengene, where MM reflects the module connectivity of each gene. GS refers to the correlation coefficient between genes and clinical traits, representing the correlation between each gene and DCM-induced HF (20). Hence, genes with larger absolute GS and MM values were associated with DCM-induced HF. Genes with |MM| >0.8 and |GS|>0.8 in the clinically relevant gene modules were defined as hub genes. Finally, the overlapping parts of hub genes in the crucial modules and DEGs were considered as key genes and visualized using a Venn diagram.

Validation of key genes and ROC curve analyses

The GSE116250 dataset was used to verify the expressions of the key genes. In addition, the ROC curves and the area under the curve (AUC) were performed to calculate the diagnostic value of key genes identified using ‘pROC’ (v 1.18.0) R package.

Immune infiltration analysis

To estimate the relative fractions of infiltrating immune cells in DCM-induced HF, the CIBERSORT algorithm was used to determine the characteristics of immune cell infiltrations (29).

Correlation analysis between diagnostic markers and immune cells

Spearman's correlation analysis was performed to examine the correlation between immune cells and diagnostic markers.

Validation of the identified key genes in vitro

AC16 is a proliferating human cardiomyocyte cell line from human ventricular tissue. They can be differentiated in vitro and used to study molecular mechanism of cardiomyocytes in physiological and pathological settings (30). The AC16 cells were purchased from the Shanghai EK-Bioscience Biotechnology Co., Ltd. The cells were cultivated in a humidified atmosphere with 5% CO2 at 37˚C. The cells were cultured in six-well plates in DMEM (MilliporeSigma) containing 8% FBS (Gibco; Thermo Fisher Scientific, Inc.) and 1% penicillin-streptomycin solution (MilliporeSigma). Doxorubicin (DOX; 2 µM; MedChemExpress) was used to stimulate AC16 cells at 37˚C in a humidified atmosphere of 5% CO2 for 24 h.

Reverse transcription-quantitative PCR (RT-qPCR) analysis

Total RNA was isolated from AC16 cells using the RNA Rapid Extraction Kit (ShangHai YiShan Biotechnology Co., Ltd). PrimeScript™ RT Master Mix (cat. no. RR036A; Takara Bio, Inc.) was used to synthesize cDNA according to the manufacturer's protocol. SYBR Premix (cat. no. RR420A; Takara Bio, Inc.) was used for qPCR (reaction conditions: 95˚C pre-denaturation for 30 sec, 95˚C denaturation for 5 sec and 60˚C annealing for 31 sec, for 40 cycles). The relative expression levels of genes were assessed using the 2-ΔΔCq method; GAPDH was used as an endogenous loading control (31). The primers used are listed in Table SI.

Western blotting

AC16 cells were cultured in 6-well plates to 80% density and then incubated with 2 µM DOX at 37˚C for 24 h. Subsequently, the total cell proteins were extracted from AC16 cells by RIPA Lysate (Beyotime Institute of Biotechnology). SDS-PAGE gel preparation kit (cat. no. P0012A; Beyotime Institute of Biotechnology) was used to prepare the gel (5% stacking gel, 12% separating gel concentration), and then 20 µg protein per lane was added for electrophoresis on 0.22-µM PVDF membranes (MilliporeSigma). The membranes were blocked with 5% skimmed milk for 2 h at room temperature, and then incubated in anti-BCL2 (cat. no. 26593; 1:2,500; Proteintech Group, Inc.), anti-BAX (cat. no. 50599; 1:1,000; Proteintech Group, Inc.), anti-atrial natriuretic peptide (ANP; cat. no. 27426; 1:1,000; Proteintech Group, Inc.) and HRP-conjugated β-actin (cat. no. HRP-60008; 1:1,000; Proteintech Group, Inc.) overnight at 4˚C. HRP-conjugated Affinipure Goat Anti-Rabbit IgG (cat. no. SA00001-2; 1:5,000; Proteintech Group, Inc.) was added at room temperature for 2 h. After washing, proteins were visualized using an ECL luminescence kit (cat. no. WBKLS0500; MilliporSigma).

TUNEL assay

AC16 cells were cultured in 24-well plates to 80% density and then incubated with 2 µM DOX for 24 h. Apoptosis was also evaluated with a fluorescence microscope (five fields of view) using One Step TUNEL Apoptosis Detection Kit (cat. no. C1086; Beyotime Institute of Biotechnology) according to the manufacturer's instructions.

Cell Counting Kit-8 (CCK-8) assay

AC16 cells were inoculated into 96-well plates at a density of 1x104 cells/well (100 µl/well). The cells were cultured with various concentrations of DOX (0, 0.5, 1, 2 and 4 µM/ml) at 37˚C for 24 h. Next, the cells were incubated with 10 µl CCK-8 reagent (Dojindo Laboratories, Inc.) at 37˚C for 1 h in a humidified CO2 incubator. The absorbance (optical density) value was analyzed at 450 nm, according to the manufacturer's instructions.

Statistical analysis

Data are presented as the mean ± SD; Microsoft Excel software 2019 (Microsoft Corporation) and GraphPad Prism 8 software (GraphPad Software; Dotmatics) were used for data analysis. The diagnostic value of key genes was assessed by ROC curve. The differences between two groups were assessed by unpaired Student's t-test, and one-way ANOVA with Bonferroni's multiple comparison post hoc test was used to compare multiple groups. P<0.05 was considered to indicate a statistically significant difference.

Results

Gene co-expression module construction using WGCNA

From the 20 samples (9 patients with DCM-induced HF and 11 control patients) in the GSE79962 dataset, the top 4,343 genes were identified and used to construct a co-expression network. The results of the cluster analysis of the samples are presented in Fig. 2A, which shows a clear distinction between the samples from patients with DCM-induced HF and the control group. An appropriate softthresholding was screened out by analysis of scale independence and mean connectivity for various soft-threshold powers (β). The soft-thresholding power (β) was set as 9 (R2=0.9) to construct a scale-free co-expression network (Fig. 2B). MergeCutHeight was used as the dendrogram cut height for module merging. Modules with similar expression patterns were merged. By setting the mergeCutHeight as 0.25, 24 modules were obtained (Fig. 2C and D).

Figure 2 Construction of gene co-expression modules. (A) Clustering dendrograms of 20 samples Red, 11 patients with dilated cardiomyopathy-induced heart failure; green, 9 healthy controls. (B) Finding appropriate soft-threshold powers (β). (C) Clustering of modules with similar expression profile. (D) Clustering dendrogram of genes by hierarchical clustering.

Identification and functional annotation of the blue module corresponding to DCM-induced HF

The aforementioned data indicated that the blue module showed the strongest correlation with DCM-induced HF (r=0.91; P<0.001; Fig. 3A and B). Therefore, 602 genes in the blue module were extracted to perform functional annotation analysis. KEGG pathway analysis revealed that that the blue module was enriched in the ‘p53 signaling pathway’, ‘MAPK signaling pathway’, ‘AGE-RAGE signaling pathway in diabetic complications’, ‘Adrenergic signaling in cardiomyocytes’, ‘JAK/STAT signaling pathway’ and ‘cGMP/PKG signaling pathway’ (Fig. 3C). In addition, GO enrichment analysis indicated that the most significant functional terms in the blue module were mainly ‘cardiac muscle tissue development’, ‘response to oxidative stress’, ‘muscle contraction’, ‘focal adhesion’ and ‘phosphoric ester hydrolase activity’ (Fig. 3D).

Figure 3 Key module selection and functional enrichment annotation. (A) Identification of the key modules associated with DCM-induced HF. Each row refers to a specific module. P-values and the corresponding correlation are stated in the cell. Blue represents negative correlation and red represents positive correlation. (B) Gene significance for DCM-induced HF in the blue module. (C) Kyoto Encyclopedia of Genes and Genomes pathway analysis of genes involved in the blue module. (D) GO functional term enrichment analysis of genes in the blue module. BP, biological process; CC, cellular component; DCM, dilated cardiomyopathy; HF, heart failure; GO, gene ontology; MF, molecular function.

GSEA

To further explore the key biological pathways underlying the development of DCM-induced HF, GSEA was performed on the GSE79962 dataset. The results showed that the pathways of ‘AGE-RAGE signaling pathway in diabetic complications’, ‘Wnt signaling pathway’, ‘Th1 and Th2 cell differentiation’ and ‘ECM-receptor interaction’ were enriched in patients with DCM-induced HF (Fig. 4).

Figure 4 Gene set enrichment analysis in the GSE79962 dataset from the Gene Expression Omnibus database. (A) AGE/RAGE signaling pathway in diabetic complications. (B) Th1 and Th2 cell differentiation. (C) Wnt signaling pathway. (D) ECM-receptor interaction. AGE, advanced glycation end product; ECM, extracellular matrix; RAGE, AGE receptor.

Key genes identification

The blue module, which contained 602 genes, was regarded as the most relevant to DCM-induced HF compared with the other modules. Using |MM|>0.8 and |GS|>0.8 as cut-off values, 40 genes were identified as hub genes (Table SII). The DEGs in patients with DCM-induced HF were then screened with 40 genes being defined as significant DEGs (Table SIII); volcano plots show the significant DEGs (Fig. 5A). Subsequently, A Venn diagram identified eight overlapping hub genes from the WGCNA analysis and the DEGs analysis, which were considered to be key genes (Fig. 5B). The mutual genes included secreted protein acidic and rich in cysteine (SPARC)-related modular calcium-binding protein 2 (SMOC2), serpin family A member 3 (SERPINA3), myosin heavy chain 6 (MYH6), S100 calcium binding protein A9 (S100A9), tubulin α3 (TUBA3)E, TUBA3D, lymphatic vessel endothelial hyaluronic acid receptor 1 (LYVE1) and phospholipase C ε1 (PLCE1). SMOC2, and PLCE1 were upregulated in patients with DCM-induced HF compared with normal healthy controls in the GSE79962 dataset, whereas SERPINA3, MYH6, S100A9, TUBA3E, TUBA3D, and LYVE1 were downregulated (Fig. 5C).

Figure 5 Identification of key genes. (A) Volcano plots of the significant DEGs between patients with DCM-induced HF and normal healthy control in the GEO dataset GSE79962. (B) Identification of common genes between hub genes identified in the WGCNA analysis and in the DEG analysis. (C) Comparison of expression levels of key genes between patients with DCM-induced HF and normal healthy controls in GSE79962. DCM, dilated cardiomyopathy; HF, heart failure; DEG, differentially expressed genes; FC, fold change; GEO, Gene Expression Omnibus; WGCNA, weighted gene co-expression network analysis.

Validation of key candidate genes

The mRNA expression levels of the eight key genes were verified in the GSE116250 dataset. Patients with DCM-induced HF showed increased expression levels of SMOC2 (Fig. 6A) and PLCE1 (Fig. 6B), whereas the expressions of SERPINA3 (Fig. 6C), MYH6 (Fig. 6D), S100A9 (Fig. 6E), LYVE1 (Fig. 6F), TUBA3D (Fig. 6G) and TUBA3E (Fig. 6H) were decreased. The results were similar to those exhibited in the GSE79962 dataset.

Figure 6 Key gene validation in the Gene Expression Omnibus dataset GSE116250. mRNA expression levels of (A) SMOC2 and (B) PLCE1 were significantly increased in patients with DCM-induced HF compared with the control group. Expression levels of (C) SERPINA3, (D) MYH6, (E) S100A9, (F) LYVE1, (G) TUBA3D and (H) TUBA3E were significantly downregulated in patients with DCM-induced HF. DCM, dilated cardiomyopathy; HF, heart failure; LYVE1, lymphatic vessel endothelial hyaluronic acid receptor 1; MYH6, myosin heavy chain 6; PLCE1, phospholipase C ε1; S100A9, S100 calcium binding protein A9; SERPINA3, serpin family A member 3; SMOC2, cysteine-related modular calcium-binding protein 2; TUBA3, tubulin α3.

ROC analysis

To evaluate the ability of key genes to serve as potential diagnostic biomarkers of DCM-induced HF, ROC curves were performed in the GSE116250 dataset. In the GSE116250 dataset, the eight key genes exhibited high predictive accuracy for diagnosing DCM-induced HF (Fig. 7). The AUC values of SMOC2, PLCE1 and SERPINA3 were >0.9, indicating that these three key genes carried the highest accuracies. The other five key genes (MYH6, S100A9, LYVE1, TUBA3D and TUBA3E) also showed high specificity with AUCs >0.8.

Figure 7 ROC curve analysis. The eight genes included (A) SMOC2, (B) PLCE1, (C) SERPINA3, (D) MYH6, (E) S100A9, (F) LYVE1, (G) TUBA3D, and (H) TUBA3E. AUC, area under the curve; LYVE1, lymphatic vessel endothelial hyaluronic acid receptor 1; MYH6, myosin heavy chain 6; PLCE1, phospholipase C ε1; ROC, receiver operating characteristic; S100A9, S100 calcium binding protein A9; SERPINA3, serpin family A member 3; SMOC2, cysteine-related modular calcium-binding protein 2; TUBA3, tubulin α3.

Immune cell infiltration analysis

Using the CIBERSORT model, 22 types of immune cell were identified between DCM-induced HF patients and controls. Fig. 8A demonstrates the immune cell infiltration difference from 11 non-failing donors and 9 patients with DCM-induced HF. The difference in immune cell infiltration showed that the fractions of naive B cells and CD4-memory-activated T cells in DCM-induced HF groups were higher compared with the control groups, whereas the infiltration of monocytes and plasma cells was lower (Fig. 8B and C).

Figure 8 Extent of immune infiltration between patients with DCM-induced HF and those in the control group. (A) The relative proportion of various immune cells. (B) The heat map of different infiltrating immune cells. (C) The difference of infiltrating immune cells between patients with DCM-induced HF (red) and the control group (green). *P<0.05, **P<0.01 vs. non-failing donors. DCM-induced HF, dilated cardiomyopathy-induced heart failure; ns, not significant; ?, P-value cannot be calculated.

Analysis of key genes and immune cells

A significant positive correlation was identified between SMOC2 and naive B cells (r=0.49, P=0.027; Fig. 9A), and negative correlation was found between SMOC2 and monocytes (r=-0.64, P=0.0029; Fig. 9B). PLCE1 was positively correlated with naive B cells (r=0.47, P=0.036; Fig. 9C). MYH6 had a positive correlation with monocytes (r=0.49, P=0.031; Fig. 9D). SERPINA3 was positively correlated with naive B cells (r=-0.48, P=0.031; Fig. 9E). S100A9 showed a positive correlation with monocytes (r=0.47, P=0.038; Fig. 9F). LYVE1 was negatively correlated with naive B cells (r=-0.47, P=0.037; Fig. 9G).

Figure 9 Correlations between key genes with immunocytes. Correlation between SMOC2 and (A) naive B cells and (B) monocytes. (C) Correlation between PLCE1 and naive B cells. (D) Correlation between MYH6 and monocytes. (E) Correlation between SERPINA3 and naive B cells. (F) Correlation between S100A9 and monocytes. (G) Correlation between LYVE1 and naive B cells. LYVE1, lymphatic vessel endothelial hyaluronic acid receptor 1; MYH6, myosin heavy chain 6; PLCE1, phospholipase C ε1; S100A9, S100 calcium binding protein A9; SERPINA3, serpin family A member 3; SMOC2, cysteine-related modular calcium-binding protein 2.

Validation of the identified key genes in vitro

The mRNA expression levels of selected key genes were examined by RT-qPCR in AC16 human cardiomyocyte cells treated with DOX, which is known to induce cardiac injury can be used as an in vitro model to mimic the mechanism of HF (32,33). DOX treatment downregulated cardiomyocyte viability in a concentration-dependent manner, with viability approaching 0.5 when treated with 2 µM (Fig. 10A). Therefore, 2 µM DOX was chosen for subsequent experiments. Following treatment of AC16 cells with 2 µM DOX, the protein expression level of BAX, a molecular marker of cell apoptosis, was increased in the DOX group compared with the control group (P<0.01; Fig. 10B). The expression levels of BCl2, an anti-apoptotic protein, were decreased in the DOX group compared with the control (P<0.01; Fig. 10B). Additionally, the apoptotic levels of AC16 cells were detected by TUNEL assay, which showed that DOX stimulation significantly increased the number of apoptotic AC16 cells to 10.37% (P<0.001; Fig. 10C). The mRNA expression levels of ANP (P<0.001) and brain natriuretic peptide (BNP; P<0.01) were also increased in the DOX groups compared with the controls (Fig. 10D and E). The protein expression levels of ANP were increased in the DOX group compared with the control (P<0.01; Fig. 10F). These results indicated that an in vitro model of DOX-induced cardiac injury was successfully constructed, which was used to validate the expression of hub genes.

Figure 10 Construction of a cardiac injury model in vitro. (A) Cell viability following DOX treatment at various concentrations. ****P<0.0001 vs. 0 µM. (B) Protein expression levels of BCL2 and BAX protein in DOX-induced cardiac injury. **P<0.01 vs. control. (C) Apoptotic levels of AC16 cells in DOX-induced cardiac injury were determined by TUNEL assay. ***P<0.001 vs. control. Scale bar, 50 µm. (D) mRNA expression levels of ANP in DOX-induced cardiac injury. ***P<0.001 vs. control. (E) mRNA expression levels of BNP in DOX-induced cardiac injury. **P<0.01 vs. control. (F) Protein expression levels of ANP in DOX-induced cardiac injury model. **P<0.01 vs. control. ANP, atrial natriuretic peptide; BNP, brain natriuretic peptide; DOX, doxorubicin.

As shown in Fig. 11, the gene expression for PLCE1, SERPINA3, MYH6, S100A9, and LYVE1 was similar to that of the microarray analysis in the GSE79962 and GSE116250 datasets. Specifically, PLCE1 (P<0.0001; Fig. 11B) was upregulated in DOX-induced cardiac injury, whereas SERPINA3 (P<0.0001; Fig. 11C), MYH6 (P<0.001; Fig. 11D), S100A9 (P<0.01; Fig. 11E) and LYVE1 (P<0.001: Fig. 11F) were downregulated compared with the control group. However, the expression of SMOC2, TUBA3D and TUBA3E were inconsistent with the results of microarray datasets. Specifically, the expression of SMOC2 was decreased in DOX-induced cardiac injury compared with the controls, but the difference was not statistically significant (Fig. 11A). The expression of TUBA3E was increased in DOX-induced cardiac injury compared with that in untreated cells (Fig. 11H), whereas there was no difference in the levels of TUBA3D (Fig. 11G).

Figure 11 Relative mRNA expression levels of key genes between control and DOX-induced cardiac injury cells. mRNA expression levels for (A) SMOC2, (B) PLCE1, (C) SERPINA3, (D) MYH6, (E) S100A9, (F) LYVE1, (G) TUBA3D, and (H) TUBA3E. **P<0.01, ***P<0.001, ****P<0.0001 vs. control. Control, untreated cells; DOX, doxorubicin; LYVE1, lymphatic vessel endothelial hyaluronic acid receptor 1; MYH6, myosin heavy chain 6; PLCE1, phospholipase C ε1; S100A9, S100 calcium binding protein A9; SERPINA3, serpin family A member 3; SMOC2, cysteine-related modular calcium-binding protein 2; TUBA3, tubulin α3.

Discussion

DCM is one of the leading causes of HF worldwide and is the most common indication for heart transplantation; it is characterized by an enlarged heart with impaired contractility (34). Individuals with the lowest ejection fractions, often accompanied by circulatory collapse, arrhythmias and thromboembolic events, have the worst prognosis (35). Despite the emergence of novel therapies, such as cardiac contractility modulation therapy and stem cell therapy, DCM-induced HF still has a high morbidity and mortality rate (36-40). Therefore, elucidating the underlying mechanism of DCM-induced HF may aid in the diagnosis and treatment of patients. The present study aimed to investigate the gene co-expression networks in DCM-induced HF. Using the WGCNA method, the blue module was selected as the key module corresponding to DCM-induced HF, and pathway enrichment analysis revealed that the AGE/RAGE signaling pathway in diabetic complications, the p53 and MAPK signaling pathways, adrenergic signaling in cardiomyocytes, and the JAK-STAT and cGMP-PKG signaling pathways were associated with DCM-induced HF.

Previous studies have shown that the AGE/RAGE signaling pathway in diabetic complications can trigger fibroblast activation in the heart, leading to fibroblast-mediated matrix remodeling in HF (41,42). The present study further emphasizes the crucial roles of the AGE/RAGE signaling pathway in diabetic complications in the pathogenesis of DCM-induced HF. Previous studies have indicated that aberrant JAK/STAT signaling may promote progression from hypertrophy to HF (43-45). In addition, the cGMP/PKG signaling pathway, adrenergic signaling in cardiomyocytes, and the p53 signaling pathway are also reportedly involved in the pathogenesis of HF (46-52). The enrichment findings of the present study are in line with the aforementioned conclusions, further strengthening the reliability of the present results.

Using WGCNA and DEGs analysis, eight key genes were identified in the blue module: SMOC2, SERPINA3, MYH6, S100A9, TUBA3E, TUBA3D, LYVE1 and PLCE1. SMOC2, and PLCE1 were upregulated in tissues from patients with DCM-induced HF compared with healthy controls, whereas the expression of SERPINA3, MYH6, S100A9, TUBA3E, TUBA3D and LYVE1 was lower. Furthermore, the expression of these key genes was validated in the GSE116250 dataset. The expression levels of SERPINA3, MYH6, S100A9, LYVE1 and PLCE1 were also validated in vitro. The present study found that these genes may perform a crucial role in the pathophysiology of DCM-induced HF.

The present study not only identified the MYH6 gene already involved in DCM and HF, but also provided new candidate genes (SERPINA3, SMOC2, S100A9, LYVE1 and PLCE1) for further experimental investigation. MYH6, encoding the α heavy chain subunit of cardiac myosin, is expressed primarily in human cardiac atria and performs crucial roles in cardiac muscle contraction (53). Mutations in the MYH6 gene are associated with dilated as well as hypertrophic phenotypes of cardiomyopathy (54-57). The subsequent structural changes of the myocardium induced by alterations in MYH6 expression may eventually lead to cardiac enlargement and dysfunction (58). SERPINA3 is a member of the serpin superfamily of protease inhibitors involved in a wide range of biological processes (59). Previous studies have found that SERPINA3 is downregulated in DCM and HF, and plasma SERPINA3 levels are associated with poor survival in patients with HF (60-63). It is mainly involved in regulating the inflammatory response and oxidative stress (64,65); however, the function of SERPINA3 in HF is unknown. We hypothesized that SERPINA3 may participate in the progression of HF by regulating inflammatory activity. SMOC2, a member of the SPARC family of matricellular proteins, modulates cell-matrix interactions (66). It has been previously reported that SMOC2 expression is higher in right ventricular failure tissue (67). In addition, high expression of SMOC2 is associated with cardiac fibrosis in chronic Chagas disease cardiomyopathy (68). A number of studies have reported that SMOC2 promotes tissue fibrosis by regulating fibroblast-to-myofibroblast transformation (66,69,70). As cardiac fibrosis is important in HF progression (71), it is reasonable to suggest that high expression of SMOC2 may contribute to cardiac fibrosis in DCM-induced HF. S100A9, a Ca2+-binding protein, has anti-inflammatory and immunoregulatory actions, it serves key immune response roles in inflammatory disorders, including cardiovascular disease (72-75). A study showed that recombinant S100A8/A9 attenuates cardiac hypertrophy and fibrosis by suppressing the calcineurin/NFAT pathway (76). Marinkovic et al (77) found that short-term S100A9 blockade reduces cardiac inflammation, limits myocardial damage, and significantly improves cardiac function and hemodynamics following myocardial ischemia; however, long-term S100A9 blockade negatively impacts cardiac recovery. Marinkovic et al (78) further found that short-term S100A9 blockade improves cardiac function following myocardial infarction by inhibiting inflammation. This suggested that the downregulation of S100A9 may be closely associated with the occurrence of DCM-induced HF. LYVE-1, a docking receptor for hyaluronic acid-coated leukocytes, regulates the activation of lymphocytes and the entry of immune cells from tissues into lymphatic vessels (79,80). A previous study in mice reported that LYVE-1 deletion leads to increased chronic inflammation and long-term deterioration of cardiac function (80). PLCE1, a member of the phosphoinositide-specific PLC family, is essential for intracellular signaling through the catalyzation of membrane phospholipid hydrolysis (81). Overexpression of PLCE1 has been reported to promote inflammation in myocardial ischemia-reperfusion injury through the activation of the NF-κB signaling pathway (82). However, the associations between these key genes and the mechanisms of DCM-induced HF have not been identified and are worth exploring in the future.

The diagnostic value of these genes was also explored in the present study. The eight key genes showed a robust predictive value in DCM-induced HF, indicating their potential use as biomarkers. In addition, the characteristics of immune cell infiltrations in DCM-induced HF were also explored, which revealed a significant enrichment of naive B cells as well as CD4-memory-activated T cells in the DCM-induced HF samples. Previous studies have shown that CD4+ T cells are abnormally activated in patients with DCM and have a direct pathogenic role in the development of HF (83-85).

The present study has some limitations, including the failure to assess BNP, NT-proBNP, TnI and TnT using western blotting for in vitro phenotype validation. Additionally, the present study did not validate the bioinformatics results in DCM-induced HF and normal human tissues in vivo. Although eight key genes associated with DCM-induced HF were identified, the specific mechanism of these genes was not demonstrated. In future, in vivo experiments will be performed to explore the specific effects of these genes in DCM-induced HF.

In conclusion, the current study identified the functional pathways, infiltrating immune cells and key genes associated with DCM-induced HF, which may aid in our understanding of the pathology and molecular mechanism of DCM-induced HF.

Supplementary Material

The primers used in this study.

Hub genes in WGCNA analysis.

Significant DEGs in GSE79962.

Acknowledgements

Not applicable.

Funding

Funding: This study was funded by The Natural Science Foundation of Shanghai (grant. no. 19JC1415703) and the Foundation of Jinshan Hospital of Fudan University (grant. no. JYQN-LC-202007).

Availability of data and materials

The datasets generated and/or analyzed during the current study are available from GEO at (https://www.ncbi.nlm.nih.gov/geo).

Authors' contributions

LZ and HG confirm the authenticity of all the raw data. LZ and HG were responsible for study conception, as well as writing the manuscript and performing data analyses. FP performed GO and KEGG analysis. JXL was responsible for the bioinformatic data collection, figure preparation and experimental data analysis. All authors have read and approved the final manuscript.

Ethics approval and consent to participate

Not applicable.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References