Identification and Validation of Dilated Cardiomyopathy-Related Genes via Bioinformatics Analysis

Introduction

Dilated cardiomyopathy (DCM) is a primary myocardial disease with an undetermined cause. DCM is characterized by enlargement of the left or right ventricle or bilateral ventricles with a hypo-systolic ventricular function with or without congestive heart failure. Ventricular or atrial arrhythmias are common. DCM is progressive and can result in death at any stage of the disease.^1–3 Because of the unclear mechanism underlying DCM, patients are not be treated in time, and the treatment approaches available are very limited.⁴ Thus, it is necessary to study and identify novel biomarkers, which are of great interest in the diagnosis and preventive treatment of DCM.

MicroRNAs (miRNAs) are approximately 18–25 nucleotide-long endogenously expressed non-coding single-stranded RNAs that specifically regulate the expression of target genes.⁵ miRNAs have been identified to be altered in cardiac disease.⁶ It has also been reported that there is a specific relationship between the dysregulation of miRNAs and cardiac development and heart failure.^7–9 Paweł Rubiś et al found that miR-133a is related to cardiovascular events (DCM).¹⁰ Furthermore, Kenta Hirai et al confirmed that cardiosphere-derived exosomal miRNAs are related to DCM.¹¹ In previous studies, it was elucidated that the 3’-UTR of miRNAs can be bound by multiple RISC complexes, which may lead to synergistic co-repression of the target mRNA. Moreover, it has been found that complexes on the 3’-UTR of miRNAs may alter the structural organization of the target mRNAs, thereby affecting the progression of various Diseases.^12,13

With the development of public databases, the gene chip microarray is an open-access database that was used to further investigate the biological process and function of various diseases. In the present study, 3 microarrays were chosen from the Gene Expression Omnibus (GEO) database to analyze the differentially expressed genes, gene-gene interaction, pathway, and GO enrichment in DCM. In addition, we also predicted several target genes and drugs for the treatment of DCM. The flow chart of the study is presented in Figure 1.

Figure 1 The flow chart of the study.

In our study, we observed differentially expressed miRNAs (DEMis) and mRNAs (DEMs) in DCM tissues and normal tissues. Therefore, we hypothesized that DEMis and DEMs may be involved in the pathogenesis of DCM. To explore the underlying mechanism of DCM, difference analysis and functional analysis were performed. Experiments were performed to verify the differential expression of hub genes. The prediction of drugs in DCM may be a good reference in future studies. Thus, the data obtained in the present study may help find novel biomarkers and therapeutic targets for DCM.

Materials and Methods

Data Resources and Processing

Microarray data from miRNA (GSE112556) and mRNA (GSE29819 and GSE84796) were obtained from the public GEO database (http://www.ncbi.nlm.nih.gov/geo/). The miRNA dataset (GSE112556) contained 3 DCM tissues and 3 normal tissues, and the expression matrix was normalized by the robust multichip average (RMA) algorithm (Supplement Figure 1). The GSE29819 dataset contained 14 DCM tissues and 12 normal tissues was analyzed with MAS 5.0 using Affymetrix default analysis settings and global scaling as the normalization method (Supplement Figure 2). The GSE84796 dataset contained 10 DCM tissues and 7 normal tissues was analyzed and normalized using GeneSpring software (11.5.1) (Supplement Figure 3).

Screening of Differentially Expressed Genes (DEG)

The “limma” R package was used to identify the differentially expressed miRNAs and mRNAs.¹⁴ P < 0.01 and |log2FC| > 1 were the screening criteria for the GSE 112556, and the adjusted P < 0.05 and |log2FC| > 1 were the screening criteria for GSE21819 and GSE84796. The differentially expressed genes volcano plots were visualized by the “ggplot2” R package and hierarchical clustering heatmaps were visualized by “Heatmap” R packages.¹⁵

Pathway and Gene Ontology Enrichment Analysis

For Pathway and GO enrichment analysis, the mRNAs were imported into the KOBAS 3.0 database (http://kobas.cbi.pku.edu.cn/). The pathway analysis terms from the Kyoto Encyclopedia of Genes and Genomes (KEGG) database were summarized by the KOBAS platform. P < 0.05 was considered statistically significant.

Construction of the miRNA-mRNA Network

The differentially expressed miRNAs were uploaded to the miRWalk v2.0 database for predicting the potential targets of miRNAs (http://mirwalk.umm.uni-heidelberg.de/).¹⁶ Then, 39 mRNAs were obtained by taking the intersection between up-regulated miRNA targets (downregulated miRNA targets) and up-regulated DEMs (downregulated DEMs). To construct the initial miRNA-mRNA network, the miRNA-mRNA pairs were imported into Cytoscape software. All mRNAs were uploaded to the STRING database (https://string-db.org/), which was used to construct a protein-protein interaction network (PPI). The initial miRNA-mRNA and PPI network were merged into the final network.

Identification of Hub Genes

The Cytohubba plugin of the Cytoscape software was used to screen the top 10 hub genes of the miRNA-mRNA network using 5 different algorithms. The 5 algorithms were degree, betweenness, closeness, stress, and Maximum Neighborhood Component (MCC), respectively. The final hub genes were obtained by taking the intersection of the results of the 5 algorithms. A Venn diagram was drawn to exhibit the overlapping hub genes of the results of the 5 algorithms.

Construction of the Cardiomyocyte Injury Model

H9c2 cells were derived from rat embryonic heart tissue and could maintain their cardiac myogenic cell characteristics. For the myocardial injury model, H9c2 cells were cultured in Dulbecco’s Modified Eagle’s Medium (DMEM) containing Doxorubicin (DOX) (2μM) (Solaribio, Life Science, Beijing, China) for 24 hours.¹⁷ The cardiomyocyte injury model was used for the following experiments (qRT-PCR).

Real-Time Quantitative Reverse Transcription PCR

Total cell RNA was extracted by using the TRIzol reagent (Invitrogen, CA, USA). Reverse transcription and the PCR process were performed according to the manufacturer’s instructions (Applied Biosystem, ThermoFisher Scientific, Shanghai, China). The mRNA primers were designed and synthesized by Sangon Biotech Co. Ltd. (Shanghai, China), and miRNA primers were synthesized by RiboBio Co. Ltd. (Guangzhou, China). The primer sequences are presented in Supplementary Table 1. The target genes were quantified relative to internal control using the 2^−ΔCT algorithm (GAPDH/U6), and the differences between the different groups were confirmed.

Statistical Analysis

The Student’s t-test was used to analyze the significant difference between the two groups. The Benjamini and Hochberg method was used to calculate the adjusted P-value for controlling the error rate. An adjusted P-value <0.05 was considered statistically significant.

Results

Because of the development of various public databases and microarray analysis, our study aimed to identify novel biomarkers. Therefore, to explore the underlying mechanisms in DCM, difference analysis and enrichment analysis were performed. Moreover, the hub gene-related drugs that were predicted may help find several novel treatments. Experiments were performed to increase the reliability of these analyses. Furthermore, our study may serve as a reference for advancing the development of DCM treatments.

Differentially Expressed miRNAs in Dilated Cardiomyopathy

The miRNA data was downloaded from the GEO database, which included 3 DCM tissues and 3 normal heart tissues. The screening criteria were described as before (|log2FC| > 1, P < 0.01). In addition, compared to normal heart tissues, 2 up-regulated (hsa-miR-770-5p and has-miR-21-5p) and 6 downregulated miRNAs (hsa-miR-144-5p, hsa-miR-9-3p, hsa-miR-144-3p, hsa-miR-451a, hsa-miR-363-3p, hsa-miR-551b-3p) were identified in the DCM tissues. The 8 miRNAs are listed in Table 1. A volcano plot was drawn to show the differentially expressed miRNAs (DEMis) (Figure 2A), and a hierarchical clustering heatmap was drawn to show the expression level of DEMis (Figure 2B).

Table 1 Differentially Expressed miRNAs (DEMis)

Figure 2 The screening and functional annotation of DEMis in DCM. (A) The 8 DEMis are shown in the volcano plot (|log2FC| > 1, P-value < 0.01). (B) The heat map shows the expression level of differentially expressed miRNA. (C) The transcription factors enrichment analysis of DEMis in DCM. (D–F) The GO functional annotation of DEMis in DCM.

Functional and GO Annotation of DEMis

The 8 differentially expressed miRNAs (DEMis) were imported into FunRich 3.1.3 software. Next, transcription factors (TFs) enrichment analysis was performed to screen several TFs that were closely related to DEMis. The results showed that the top 10 TFs, including EGR1, SP4, SP1 HOXD8, NKX6-1, HOXA5, POU2F1, IRF1, HOXA9, and HOXB4 may have a regulatory relationship with DEMis (Figure 2C). Furthermore, GO enrichment analysis was performed to identify the biological functions of DEMis. The GO biological process terms indicated that miRNAs were mainly involved in the regulation of nucleobase, nucleoside, nucleotide and nucleic acid metabolism (21.6%) (Figure 2D). The GO cellular component terms showed that most miRNAs were located in the nucleus (52.7%) and the cytoplasm (45.7%) (Figure 2E). Furthermore, the GO Molecular Function terms indicated that miRNAs were mainly enriched in transcription factor activity (10.8%) (Figure 2F).

Differentially Expressed mRNAs in Dilated Cardiomyopathy

In this study, we chose two microarrays associated with DCM from the GEO database (GSE29819 and GSE84796). Then, the “limma” R package was used to analyze the different gene expressions (Supplement Tables 2 and 3). The screening criteria were as follows: |log2FC| > 1 and adjusted P < 0.05. As shown in Figure 3A, the volcano plot indicated that 209 up-regulated and 265 downregulated mRNAs were identified as differentially expressed RNAs in GSE29819. In addition, the volcano plot in Figure 3C was drawn to show the differentially expressed RNAs in GSE84796 (including 1059 up-regulated and 398 downregulated mRNAs). Moreover, the hierarchical clustering heatmap was drawn to show the expression level of DEMis (Figure 3B and D). Finally, Venn diagrams were drawn to show the overlapping up-/down-regulated mRNAs (49 up-regulated and 43 downregulated mRNAs) between GSE29819 and GSE84796 (Figure 3E and F).

Figure 3 The screening of DEMs in DCM. (A) The volcano plot was drawn to show the 209 up-regulated and 265 downregulated mRNAs in GSE29819 (|log2FC| > 1, adjusted P-value < 0.05). (B) The heatmap shows the expression level of DEMs in GSE29819. (C) The volcano plot was drawn to show the 1059 up-regulated and 398 downregulated mRNAs in GSE84796 (|log2FC| > 1, adjusted P-value < 0.05). (D) The heatmap shows the expression level of DEMs in GSE84796. (E and F) The Venn diagrams were drawn to show the overlapping up-/down-regulated mRNAs between GSE29819 and GSE84796.

KEGG Pathway and GO Enrichment Analysis of DEMs

For investigating the underlying mechanism and function of DEMs, we uploaded 92 mRNAs to the KOBAS 3.0 database. For KEGG pathway enrichment analysis, we found that DEMs were mainly enriched in the neuroactive ligand-receptor interaction, estrogen signaling pathway, cytokine-cytokine receptor interaction, Jak-STAT signaling pathway, and pathways related to cancer, etc. A bubble plot was drawn to show the top 10 pathways from the results of KEGG analysis (Figure 4A). Using GO functional and enrichment analysis, DEMs were mainly located in the plasma membrane (GO:0005886), extracellular region (GO:0005576), cytoplasm (GO:0005737), extracellular space (GO:0005615), membrane (GO:0016020), and extracellular vesicles (GO:1903561). Besides, the biological process and molecular function of DEMs were mainly involved in protein binding (GO:0005515), collagen-containing extracellular matrix (GO:0062023), identical protein binding (GO:0042802), and G protein-coupled receptor signaling pathway (GO:0007186) (Figure 4B).

Figure 4 Pathway and GO enrichment analysis of DEMs. (A) The top 10 pathways from the results of KEGG analysis (P < 0.05). (B) The top 10 terms from the results of GO enrichment analysis (P < 0.05).

Construction of miRNA-mRNA Network and Identification of Hub Genes

In this study, we first identified 39 genes by taking the intersection between up-regulated miRNA targets (downregulated miRNA targets) and up-regulated DEMs (downregulated DEMs) (Supplement Figure 4). Next, the DEMis and DEMs were imported to construct the initial network. Moreover, DEMs were uploaded to the STRING database for constructing a protein-protein interaction network (PPI). The PPI network was also imported into the Cytoscape (version: 3.8.2) software. Subsequently, the initial network and PPI network were merged into the final miRNA-mRNA network using the Cytoscape software (Figure 5A). Hub genes were identified by using the Cytohubba plugin of the Cytoscape software. Next, we chose 5 algorithms to screen the hub genes of the PPI network. As shown in Figure 5B–F, the top 10 genes were identified by degree, MCC, betweenness, closeness, and stress. Lastly, 6 genes (hsa-miR-551b-3p, hsa-miR-770-5p, PIK3R1, DDIT4, CXCR4, and hsa-miR-363-3p) were identified as hub genes by taking the intersection of 5 algorithms in the Venn diagram (Figure 5G) (Table 2).

Table 2 The Five Algorithmic Scores of Hub Genes

Figure 5 Construction of miRNA-mRNA network. (A) The network of DEMis and DEMs. (B) The top 10 hub genes from the results of Betweenness algorithms. (C) The top 10 hub genes from the results of Closeness algorithms. (D) The top 10 hub genes from the results of Degree algorithms. (E) The top 10 hub genes from the results of MCC algorithms. (F) The top 10 hub genes from the results of Stress algorithms. (G) The Venn diagram was drawn to show the overlapping hub genes in the results of 5 algorithms.

Prediction of Candidate Drugs

The DSigDB database was used to predict several drugs which are related to the 3 DEMs. The top 10 candidate drugs related to the DEMs were chosen according to the P-value and adjusted P-value. Table 3 shows the top 10 candidate drugs from the DSigDB database and the 10 drugs are shown in Figure 6A. Figure 6B shows the structure of several drugs related to the hub gene.

Table 3 Suggested Top 10 Drugs Compounds for the DCM

Figure 6 Prediction of top 10 candidate drugs. (A) The top 10 hub gene-related drugs. (B) The diagram of chemical structure of drugs.

Validation of Abnormal Hub Genes

To verify our results, we constructed a DOX-induced myocardial injury model of DCM and determined the expression of differential miRNAs using real-time PCR. The results are shown in Figure 7. Compared with the control group, the expression of hsa-miR-363-3p and hsa-miR-551b-5p was down-regulated in DOX-treated H9c2 cells; however, the expression of hsa-miR-770-5p was significantly up-regulated (Figure 7A). In addition, the expression of Hub genes CXCR4, DDIT4 and PIK3R1 were significantly upregulated in H9c2 cells treated with DOX compared to the control group (Figure 7B). These findings were consistent with our predicted results.

Figure 7 Validation of abnormal hub genes. (A) The expression level of miRNAs (* P < 0.05, ** P < 0.01). (B) The expression level of mRNAs (* P < 0.05, ** P <0.01, **** P < 0.0001).

Discussion

With the development of various public databases, microarray analysis has gradually become an optimal approach to identify differentially expressed genes that may affect the progression of various diseases.^18–20 In recent years, there are an increasing number of studies related to gene chip analysis. This is because gene chip analysis not only helps to identify new biomarkers but also helps to predict treatments and prognostic factors in DCM.²¹ To date, relatively few studies have been performed on DCM. By the time DCM is diagnosed, the disease is often at a stage when surgical treatment is necessary, so it is essential to identify predictive and diagnostic biomarkers for early diagnosis.

It has been reported that miRNAs are involved in various diseases, such as cancer,²² cardiovascular diseases,^23–25 and other diseases.^26–28 We identified 3 miRNAs and 3 mRNAs as hub genes in DCM. Literature retrieved from PubMed showed several studies about the genes we identified before, but the underlying mechanisms of action and functions in DCM remain unclear. In previous studies, it has been shown that hsa-miR-551b may regulate the progression of lung adenocarcinoma and colorectal cancer.^29,30 Feng et al found that hsa-miR-551b-5p was involved in diabetic cardiomyopathy.³¹ Besides, the serum miR-551b-3p as a potential tumor suppressor was found to be downregulated in gastric cancer tissues, and it may predict the outcome of gastric cancer.^32,33 Also, Chang et al revealed that miR-551b-3p inhibits the growth of cholangiocarcinoma cells via targeting Cyclin D1.³⁴ What’s more, Zhou et al reported that hsa-miR-363-3p can reduce endothelial cell inflammatory responses in coronary heart disease.³⁵ And has-miR-363-3p as an tumor suppressor affects the progression of multiple cancers, including osteosarcoma,³⁶ hepatocellular carcinoma,³⁷ papillary thyroid carcinoma,³⁸ and so on. All these studies suggest that has-miR-363-3p and hsa-miR-551b-3p have the function of inhibiting proliferation and affecting cell cycle. Therefore, these results provide some theoretical support for our study. And it reveals that the dysregulation of the two miRNAs may affect the iterative renewal of cells, leading to a bad phenotype of H9c2 cells. Then, Satake et al found that hsa-miR-770-5p was associated with hyperglycemia in patients with type 1 diabetes.³⁹ Lee et al also found that hsa-miR-770-5p can sensitize breast cancer to radiation.⁴⁰ And these two reports show that hsa-miR-770-5p may induce the death of cells. Therefore, studying miRNAs in DCM disease is feasible.

In addition to miRNAs, the hub mRNAs that we identified before have been studied in other diseases, which can be a reference for our research. The C-X-C motif chemokine receptor 4 (CXCR4) encodes a CXC chemokine receptor that is specific for stromal cell-derived factor-1 and can affect disease progression. Chu et al found that CXCR4 antagonism can reduce cardiac fibrosis and improve cardiac performance in DCM.⁴¹ Besides, Damås et al elucidated that CXCR4 may enhance the progression of heart failure.⁴² More importantly, Sara Beji found that CXCR4 expression is upregulated in DOX-induced cardiomyocyte injury.⁴³ DNA damage-inducible transcript 4 (DDIT4) was confirmed to promote gastric cancer and acute myeloid leukemia progression.^44,45 Previous system analysis suggested that suppression of DDIT4 can preserve heart contraction.⁴⁶ Furthermore, DDIT4 also can control cell death in cardiac I/R injury by regulating several apoptosis-related signaling pathway.⁴⁷ In addition to cardiac diseases, relevant studies in other diseases have shown that DDIT4 mediates the apoptotic process, for example, Lynn M Knowles et al verified that increasing DDIT4 induces caspase-8 mediated tumor cell apoptosis via negatively regulating the mTOR signaling pathway.⁴⁸ And Li revealed that DDIT4 can promote the production of ROS in liver injury,⁴⁹ which suggests that the expression of DDIT4 is related to the injury of cells. The detailed functional mechanisms of DDIT4 in DCM remain unclear. DDIT4 as a novel mediator of cardiomyocyte apoptosis can be a good reference for our study and future studies. Phosphoinositide-3-kinase regulatory subunit 1 (PIK3R1) was identified as an oncogene in gastric cancer, breast cancer, and pancreatic cancer.^50–52 Furthermore, PIK3R1 was found to be involved in the pathogenesis of coronary artery disease.⁵³ Yan et al have demonstrated that PIK3R1 can induce the suppression of invasiveness in breast cancer.⁵⁴ This means that the effect of PIK3R1 on the cells is controversial. And studies focused on PIK3R1 and DCM are still limited. However, Huang et al found that PIK3R1, which is up-regulated in DCM tissues, may be involved in the pathogenesis of DCM.⁵⁵ Therefore, the study of PIK3R1 in DCM may be an interesting and meaningful direction.

In the present study, we utilized bioinformatics analysis to predict potential biomarkers in DCM. Microarray analysis was used to screen DEMis and DEMs. Then, a miRNA- mRNA network was constructed, and hub genes were identified by using 5 algorithms. KEGG and GO enrichment analyses were performed to investigate the function and underlying mechanism of action of differentially expressed genes (DEGs). Moreover, the DSigDB database was used to predict several drugs, which are related to the DEGs.

Our study has some shortcomings. Due to the difficulty of obtaining samples in the short term, we only conducted simple cell-experimental validation, but we intend to conduct a further study at a later stage when we have sufficient funds and samples. Despite the shortcomings, the data obtained is of great significance for subsequent research.

Conclusion

In summary, we comprehensively obtained several differentially expressed genes (hsa-miR-551b-5p, hsa-miR-770-5p, hsa-miR-363-3p, CXCR4, DDIT4, and PIK3R1) and several drugs that may be associated with the pathogenesis of DCM.