Bioinformatics Analysis for Identifying Micro-RNAs, Long Noncoding RNAs, Transcription Factors, Immune Genes Regulatory Networks, and Potential Therapeutic Agents for Diabetic Cardiomyopathy Using an Integrated Bioinformatics Analysis


 Objectives

We identified functional genes and studied the underlying molecular mechanisms of diabetic cardiomyopathy (DCM) using bioinformatics tools.
Methods

Original gene expression profiles were obtained from the GSE21610 and GSE112556 datasets. We used GEO2R to screen the differentially expressed genes (DEGs). Gene Ontology and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analyses were performed on DEGs. Protein–protein interaction (PPI) networks of DEGs were constructed using STRING and hub genes of signaling pathways were identified using Cytoscape. Aberrant hub gene expression was verified using The Cancer Genome Atlas dataset. Connectivity Map was used to predict the drugs that could treat DCM.
Results

The DEGs in DCM were mainly enriched in the nuclei and cytoplasm and involved in DCM- and chemokine-related signaling pathways. In the PPI network, 32 nodes were chosen as hub nodes and an RNA interaction network was constructed with 517 interactions. The expression of key genes (JPIK3R1, CCR9, XIST, WDFY3.AS2, hsa-miR-144-5p, and hsa-miR-146b-5p) was significantly different between DCM and normal tissues. Danazol, ikarugamycin, and semustine were identified as therapeutic agents against DCM using CMAP.
Conclusion

The identified hub genes could be associated with DCM pathogenesis and the above drugs could be used for treating DCM.


Introduction
Diabetic cardiomyopathy (DCM) is a speci c cardiac condition in patients with diabetes, and is caused by hyperglycemia and insulin resistance. DCM development is independent of coronary artery disease, hypertension and other usual cardiac risk factors. It is clinically characterized by left ventricular hypertrophy, diastolic function damage, cardiomyocyte hypertrophy, myocardial brosis, and cardiomyocyte apoptosis, eventually leading to heart failure 1 . The pathophysiological mechanism of DCM could be related to molecular conditions such as metabolic disorders, calcium imbalance, mitochondrial damage, oxidative stress, increased in ammation, vascular dysfunction, cardiac autonomic neuropathy, micro-RNAs (miRNAs) and long noncoding RNAs (lncRNAs) dysregulation 2, 3 .
Due to the lack of clear diagnostic criteria, the speci c pathogenic mechanisms of DCM remain unclear and speci c histological characteristics, clinical manifestations and biomarkers for the de nitive diagnosis of DCM still need more investigation [4][5][6] . Therefore, screening the potential molecules and pathways implicated with DCM and further studying the molecular mechanisms of DCM is important for preventing it and developing effective treatments.
Bioinformatic analysis and gene expression pro ling analysis showed endogenous noncoding RNAs, miRNAs are widely present in eukaryotes and play an important role in the development of cardiovascular diseases 7 . miRNAs participate in cell proliferation and differentiation, metastasis, apoptosis, immune response, and other biological processes. They also play an important role in the development of DCM; more than 300 miRNAs have been con rmed to be expressed in the cardiomyocytes. Abnormal expression of miRNAs has been con rmed to be closely related to lncRNAs and transcription factors (TFs). For instance, lncRNAs can regulate miRNA expression by binding to target miRNAs and participating in the regulation of mRNA expression 8, 9 . TFs play an important role in gene transcription and post-transcriptional regulation and participate in the control of miRNA signaling pathways 10,11 .
However, further research is required to determine the regulatory mechanisms of miRNAs, lncRNAs, TFs, and mRNAs in DCM. Microarray analysis can rapidly identify all genes expressed at a given time; thus, future research can bene t from the integration and analysis of microarray data 12,13 . In this study, we aimed to conduct an in-depth investigation of the miRNAs, lncRNAs, target genes, and pathways involved in the development of DCM.

Obtaining raw data
The datasets used in the present study were downloaded from the National Center of Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) 14 . The GSE112556 dataset includes microRNA expression pro les. The original gene expression pro les were obtained from the GSE21610 and GSE17800 datasets. The GSE21610 dataset includes the gene expression pro le of DCM tissues from 42 patients with DCM and normal tissues from 8 patients (Table 1). GPL570 Affymetrix HumanGenome U133 Plus 2.0 Array was used to analyze these data. Theplatform used for these data is also speci ed in Table 1. Prediction of mRNA-miRNA-lncRNA interactions We used miRWalk 3.0 (http://mirwalk.umm.uni-heidelberg.de/) to predict the interaction between differentially expressed miRNAs and mRNAs, which integrated the prediction results of TargetScan and miRDB 16, 17 . A score ≥0.95 was considered as the cut-off criterion for predictive analysis in miRWalk. We selectedonlythe target mRNA identi ed in all the above databases for further analysis. DIANA-LncBase v2.0 18 was used to predict the interaction between miRNA and lncRNA, and a score of ≥0.4 was considered as the cut-off criterion for predictive analysis in the LncBase experimental module. After determining the number of predicted DEGs, the miRNAs, lncRNAs, and mRNAs were selected for further analysis. The Cytoscape software (version 3.40) was used to visualize the supervision network.

Gene function analysis
The DAVID database (http://david.abcc.ncifcrf.gov/) was used to perform gene ontology (GO)and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of signi cantDEGs and target genes of miRNAs with differential expression. The species was limited to Homo sapiens, and the "adjusted P value (from the Benjamini-Hochberg method), 0.05" was considered statistically signi cant 19 . The GO term includedthe following three criteria: molecular function (MF), cell composition (CC), and biological process (BP).

Protein-protein interaction (PPI) network construction and analysis
The STRING online database and PPI pairs with a combined score of ≥0.4 were used to construct a PPI network. TFs were annotated with TF checkpoints 20 . Then, the Cytoscape software (version 3.4.0) was used to visualize the regulatory relationship among genes and analyze the topological properties of the network, including the network degree distribution using the CentiScaPe application 21 . In addition, according to previous studies, genes with a degreeof>5 were de ned as pivotal genes in the regulatory network 22 . Proteins with high degrees tend to be key factors.
TCGA dataset analysis TCGA is a platform for researchers to download and evaluate free public datasets (https://cancergenome.nih.gov/) 23 . In this study, we used the online tool UALCAN (http://ualcan.path.uab.edu/) 24 to verify the expression of the hub genes in the TCGA dataset in order to improve the reliability of our analysis. UALCAN uses TCGA Level 3 RNA-seq and clinical data. As recommended by Li and Dewey 25 , the Practical Extraction and Reporting Languageprogramwas used to multiply "scaledestimate" by 10 6 to obtain the transcripts per million (TPM) expression value. Since TPM is considered more comparable than reads per kilobase transcripts per million mapped reads (FPKM) and reads per kilobase transcript reads per million mapped reads (RPKM) 26 , it was used as a measurement of gene expression in our study. Box and whisker plots were generated to represent the hub gene expression between normal samples and DCM samples in pipeline for RNA sequencing data analysis.

CMAP analysis
We screened the effective drugs used for treating DCM by using the CMAP database. The DEG le was uploaded onto the CMAP website. The screening criteria were set to a mean valueof ≤0.4 and a P value of<0.05.

Results
Identi cation of differentially expressed miRNAs, lncRNAs, TFs, and immune gene mRNAs  Table 2. Considering the BP criterion, the DEGs were mainly enriched in BPs related to the regulation of biological and cellular processes. Considering the MF criterion, the DEGs were mainly enriched in nucleic acid binding, transcription regulator activity, and signaling receptor activity. Considering the CC criterion, the DEGs were mainly enriched in cell, membrane-bounded organelle, and intracellular membrane-bounded organelle. The top 20 most signi cantKEGG pathway terms are shown in Fig. 4D. The genes were mainly enriched in the cytokine-cytokine receptor interaction, PI3K-Akt signaling pathway, and chemokine signaling pathway.

PPI analysis
To further explore the most signi cant clusters of DEGs,we performed a PPI network analysis of DimGene and DTFGene. A con dence score of ≥0.9 was set as the threshold and protein nodes that did not interact with other proteins were removed. A total of 162 of 889 DEGs were mapped into the PPI network complex and the PPI network data le of STRING was imported into Cytoscape. In this network, with a degree of >5, 32 nodes were chosen as hub nodes, including 3 TFs and 22 immune genes mRNAs. The results are presented in Fig. 5A and B. In the PPI networks, we identi ed the following top 5 DEGs with the highest node degrees: phosphoinositide-3-kinase regulatory subunit 1 (PIK3R1; degree, 24), pro-melanin concentrating hormone (PMCH; degree, 19), proenkephalin (degree, 19), pro-platelet basic protein (PPBP; degree,16), and epidermal growth factor receptor (EGFR; degree, 16). The most signi cant hub genes whose protein levels were upregulated were PIK3R1, JAK2, EGFR, IL7R, CD2, CXCR4, KLRD1, IGLL5, DAPP1, HTR2B, GNRHR, and CYSLTR1, and the most signi cant hub genes whose protein levels were downregulated were ARTN, IL5RA, NGFR, LRRFIP1, ESR2, LTB4R, GCG, PMCH, GHSR, PTPN1, BCAR1, CDH1, and IL1RAP (Fig. 5C). These ndings suggest that these nodes may play important roles in the development of DCM.
mRNA-miRNA-lncRNA network analysis The miRNA-mRNA regulatory network was constructed, and interaction analysis indicated that 13 differentially expressed miRNAs targeted 26 mRNAs. In detail, the miRNAs downregulated 10 mRNAs but upregulated 16 mRNAs. We also constructed the miRNA-lncRNA regulatory network. Interaction analysis indicated that differentially expressed miRNAs targeted 13 lncRNAs, of which 6 were upregulated and 7 were downregulated. Based on the ceRNA theory, we used the sharedmiRNA as a junction and analyzed the lncRNA-miRNA-mRNA ceRNA regulatory relationship (Table 3). In the miRNA-mRNA and miRNA-lncRNA interaction pairs, the downregulated miRNAs accompanied the upregulated lncRNAs and mRNAs, and upregulated miRNAs accompanied the downregulated lncRNAs and mRNAs. Finally, the lncRNA-miRNA-TF-immune gene pathway ceRNA network was constructed using 517 interactions, including 10 miRNAs, 2 lncRNAs, and 161 mRNAs (123 immune genes and 38 TFs) (Fig. 6A). The lncRNA-miRNA-Hub gene pathway ceRNA network of key immune genes was constructed through CytoHubba and MCODE, and two important pathways, PIK3R1 and CCR9, were obtained as shown in Fig. 6B and C. including PIK3R1, CCR9, WDFY3.AS2, XIST, hsa-miR-144-5p, andhsa-miR-146b-5p, was signi cantly different between DCM and normal tissue samples (Fig. 7).

Discussion
DCM is a serious complication associated with diabetesand increasing incidence rate; however, the molecular mechanism underlying the development of DCM is unclear. Therefore, studying its mechanism and determining molecular targets for diagnosis and treatment are essential. Some studies have found that mRNAs, miRNAs, lncRNAs, and TFs play important regulatory roles in the development of DCM 2, 3 . In this study, we rstly performed a comprehensive bioinformatics analysis to identify the miRNAs, lncRNAs, immune genes, and TFs in the interaction network and reported the key genes related to DCM.
Noncoding RNAs (ncRNAs) have been proved to be of great signi cance in human life, health, and disease diagnosis. According to their functions, ncRNAs are classi ed as miRNAs, lncRNAs, and circular RNAs (circRNAs) 28 . lncRNAs act as miRNA sponges that can regulate miRNA abundance and compete with mRNA for miRNA binding. miRNAs are considered to be potential tissue-speci c biomarkers and play an important role in the pathogenesis of myocardial remodeling 29 . They also contribute to the development of DCM by regulating oxidative stress, in ammation, cardiomyocyte apoptosis, and mitochondrial function 2 . Previous study has reported that miRNA-195 inhibitors can reduce diabetic myocardial hypertrophy and improve cardiac function by reducing oxidative damage, inhibiting cell apoptosis, and promoting angiogenesis 30 . Moreover, miRNAs are involved in the regulation of multiple biological functions in the mitochondria by binding to their target genes 31 . The miR-30 family, including miR-30a, miR-30b, miR-30c, and miR-30d, is highly expressed in cardiomyocytes 32 and inhibits the expression of P53 to prevent mitochondrial division and cell apoptosis. The mechanisms by which miRNAs regulate the expression of genes involved in myocardial diseases and their effect on the development of diabetes have not yet been clearlyunderstood.
In this study, we analyzed the DEGs identi ed from the GSE21610 and GSE112556 datasets, which consist TCGA dataset analyses indicated that aberrant expression of the hubgenes, including WDFY3.AS2, XIST, hsa-miR-144-5p, and hsa-miR-146b-5p, was signi cantly different between DCM and normal tissuesamples. Previous studies have shown that MiR-144 can target Nrf2 directly and down-regulated miR-144 is found to regulate oxidative stress in diabetic cardiomyocytes 33 . MiR-146a plays a regulatory role with the NF-κB signaling pathway components, which is the key mediator of in ammation and hyperglycemia 34,35 . This is consistent with our nding. The roles of different miRNAs and lncRNAs, especially WDFY3.AS2, XIST, hsa-miR-144-5p, and hsa-miR-146b-5p, in the context of DCM should be further studied, including type 1 and type 2 diabetes, the contribution of pathophysiological mechanisms including in ammation, apoptosis, hypertrophy and brosis, and oxidative stress to the development of DCM.
Our dataset analyses also indicated the aberrant expression of the hub genes PIK3R1 and CCR9 between DCM and normal tissue. The enrichment results of PIK3R1 indicated that it was signi cantly enriched in the JAK-STAT signaling pathway, T Toll-like receptor signaling pathway, TNF signaling pathway, uid shear stress, and atherosclerosis AGE-RAGE signaling pathway (which cause diabetic complications), and the chemokine signaling pathway. The JAK/STAT signaling pathway is an intracellular signaling pathway closely related to cardiac hypertrophy, which plays a key role in cell growth, survival and differentiation, and regulation of gene expression. It has been shown that high glucose activates the JAK/STAT signaling of vascular endothelial cells in vitro through phosphorylation of JAK2 and subsequently STAT3, leading to the proliferation of endothelial cells. Hyperglycemia, hyperlipidemia, and insulin resistance are the main causes of chronic low-grade in ammation of the diabetic heart. Cardiac Toll-like receptors and in ammasome complexes, probably through NF-B activation and ROS overproduction, may be key inducers for in ammation 36 . Therefore, what role of PIK3R1 plays in regulation of JAK/STAT signaling pathway, T Toll-like receptor signaling pathway and others represents an important strategy for the treatment of DCM and needs to be further explored.
We also performed CMAP analysis to examine the potential small-molecule drugs that can be potentially used for treating DCM. The three most signi cant drugs identi ed were danazol, ikarugamycin, and semustine. Danazol is a medicine used to treat endometriosis, brocystic breast disease, hereditary angioedema and other diseases. Ikarugamycin is a previously discovered antibiotic which has been shown to inhibit clathrin-mediated endocytosis and the uptake of oxidized low-density lipoprotein in macrophages 37 . Semustine is a drug used in chemotherapy 38 . Although there is no relevant literature regarding the effect of these compounds on DCM treatment, we speculate that these small-molecule drugs may have the potential for effectively treating DCM by interfering with the expression of the hub genes.

Conclusion
We constructed and analyzed mRNAs, miRNAs, lncRNAs, and TF interaction networks to identify the key genes related to the development of DCM. Our ndings suggest that DEGs were mainly involved in the cytokine-cytokine receptor interaction, chemokine signaling pathway, and PI3K-Akt signaling pathway, and the hub genes JPIK3R1, CCR9, XIST, WDFY3.AS2, hsa-miR-144-5p, and hsa-miR-146b-5p may be the key factors associated with the pathogenesis of DCM. Danazol, ikarugamycin, and semustine could be potential targets for performing the follow-up study of the molecular mechanism underlying the development of DCM and the therapeutic drugs that can effectively treat this condition. However, further studies are required to con rm the functions of the DEGs we found and the effectiveness of these smallmolecule drugs against DCM.

Declarations
Disclosure statement The author did not declare any con icts of interest with the current research.

Funding statement
This research was supported by the Basic Research Projects (Y20211015) organized by the Wenzhou Municipal Science and Technology Bureau.

Declaration of Competing Interest
The authors declare they have no con ict of interest.