Exploring Potential Biomarkers and Molecular Mechanisms of Ischemic Cardiomyopathy and COVID-19 Comorbidity Based on Bioinformatics and Systems Biology

Cardiovascular complications combined with COVID-19 (SARS-CoV-2) lead to a poor prognosis in patients. The common pathogenesis of ischemic cardiomyopathy (ICM) and COVID-19 is still unclear. Here, we explored potential molecular mechanisms and biomarkers for ICM and COVID-19. Common differentially expressed genes (DEGs) of ICM (GSE5406) and COVID-19 (GSE164805) were identified using GEO2R. We performed enrichment and protein–protein interaction analyses and screened key genes. To confirm the diagnostic performance for these hub genes, we used external datasets (GSE116250 and GSE211979) and plotted ROC curves. Transcription factor and microRNA regulatory networks were constructed for the validated hub genes. Finally, drug prediction and molecular docking validation were performed using cMAP. We identified 81 common DEGs, many of which were enriched in terms of their relation to angiogenesis. Three DEGs were identified as key hub genes (HSP90AA1, HSPA9, and SRSF1) in the protein–protein interaction analysis. These hub genes had high diagnostic performance in the four datasets (AUC > 0.7). Mir-16-5p and KLF9 transcription factor co-regulated these hub genes. The drugs vindesine and ON-01910 showed good binding performance to the hub genes. We identified HSP90AA1, HSPA9, and SRSF1 as markers for the co-pathogenesis of ICM and COVID-19, and showed that co-pathogenesis of ICM and COVID-19 may be related to angiogenesis. Vindesine and ON-01910 were predicted as potential therapeutic agents. Our findings will contribute to a deeper understanding of the comorbidity of ICM with COVID-19.


Introduction
COVID-19 (SARS-CoV-2) is currently wreaking havoc around the world. Lung disease, the most common complication of COVID-19, can present as a progression from an asymptomatic/mild infection to multi-organ failure leading to death [1][2][3]. Besides lung disease, cardiovascular complications are also clinically recognized as risk factors for COVID-19 [4,5]. Approximately 20-30% of patients with COVID-19 who show cardiovascular manifestations have poorer outcomes [6,7]. In China, a study of 44,672 COVID-19 cases found that the death rate from combined cardiovascular complications was 10.5% compared with an overall rate of 2.4% [8]. Cardiovascular complications can be triggered by myocardial ischemia or non-ischemic processes (e.g., myocarditis) [9], and may be related to viral damage and immune responses [10][11][12].
Ischemic cardiomyopathy (ICM) reduces or stops blood flow to the heart, leading to myocardial damage [13]. In a systematic review, Omidi et al. [14] found that cardiomyopathy was a very common cardiovascular complication of COVID-19. The autopsy pathology of patients with COVID-19 found microscopic evidence that suggested acute ischemic Ischemic cardiomyopathy (ICM) reduces or stops blood flow to the heart, leading to myocardial damage [13]. In a systematic review, Omidi et al. [14] found that cardiomyopathy was a very common cardiovascular complication of COVID-19. The autopsy pathology of patients with COVID-19 found microscopic evidence that suggested acute ischemic injury [15]. Impaired cardiac systolic function was detected in patients with COVID-19 and ICM by cardiac computed tomography [16]. However, the common pathogenesis of ICM and COVID-19 remains to be fully elucidated.
Disease-disease relationships play important roles in the pathobiological manifestations of diseases and in their precise treatment [17]. Exploring disease associations can enhance the understanding of connections between diseases, which can help in developing diagnostic, prognostic, and treatment strategies [18]. Mutations in functionally relevant genes are known to be responsible for overlapping clinical phenotypes of different diseases. Therefore, in this study, we focused on key genes in ICM that are connected to the pathogenesis of COVID-19. The workflow of this study is shown in Figure 1. We analyzed two high-throughput sequencing datasets (GSE164805 and GSE5406) that we downloaded from the Gene Expression Omnibus (GEO). We used an integrated bioinformatics and systems biology approach to identify differentially expressed genes (DEGs) and explored their molecular mechanisms in COVID-19 and ICM. These genes provide an entry point to explore the common clinical phenotypes of COVID-19 and ICM. The proteinprotein interaction (PPI) analysis identified key hub genes. To explore the diagnostic performance for these hub genes, we downloaded two other GEO datasets (GSE116250 and GSE211979) and plotted receiver operating characteristic (ROC) curves. Transcription factor (TF) and microRNA (miRNA) regulatory networks for these genes were constructed. We identified potential therapeutic agents by drug prediction, and validated them by molecular docking computer simulations.

Collection of Common DEGs for COVID-19 and ICM
In the GSE5406 dataset, 204 genes were up-regulated and 139 genes were downregulated ( Figure 2A). In the GSE164805 dataset, 9268 genes were up-regulated and 9011 genes were down-regulated ( Figure 2B). We analyzed the intersection of the DEGs from these two datasets and identified 81 common DEGs; 28 were co-up-regulated and 53 were co-down-regulated for COVID-19 and ICM ( Figure 2C,D).

Collection of Common DEGs for COVID-19 and ICM
In the GSE5406 dataset, 204 genes were up-regulated and 139 genes were down-re ulated (Figure 2A). In the GSE164805 dataset, 9268 genes were up-regulated and 90 genes were down-regulated ( Figure 2B). We analyzed the intersection of the DEGs fro these two datasets and identified 81 common DEGs; 28 were co-up-regulated and 53 we co-down-regulated for COVID-19 and ICM ( Figure 2C,D).

Enrichment Analysis of Common DEGs for COVID-19 and ICM
We performed a functional enrichment analysis of the 81 common DEGs using t DAVID database. The DEGs were enriched mainly with the gene ontology (GO) term angiogenesis under the biological process category, various extracellular components u der the cellular component category, and transcriptional regulation of endoplasmic ret ulum stress under the molecular function category ( Figure 3A-C). The DEGs were e riched mainly with Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways relat to mineral absorption and protein processing in endoplasmic reticulum. ( Figure 3D These results strongly suggest that common biological processes and pathways are i volved in the development of ICM and COVID-19.

Enrichment Analysis of Common DEGs for COVID-19 and ICM
We performed a functional enrichment analysis of the 81 common DEGs using the DAVID database. The DEGs were enriched mainly with the gene ontology (GO) terms angiogenesis under the biological process category, various extracellular components under the cellular component category, and transcriptional regulation of endoplasmic reticulum stress under the molecular function category ( Figure 3A-C). The DEGs were enriched mainly with Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways related to mineral absorption and protein processing in endoplasmic reticulum. ( Figure 3D). These results strongly suggest that common biological processes and pathways are involved in the development of ICM and COVID-19.

Construction of PPI Networks to Identify Hub Genes
Individual nodes were hidden to obtain the PPI network. The network was optimized using Cytoscape and 62 nodes were detected ( Figure 4A). The top 10 hub genes in the PPI were screened using four methods (MCC, EPC, Degree, BottleNeck) in the cytoHubba plugin ( Table 1). The intersection of the results identified three pivotal genes, HSP90AA1, HSPA9, and SRSF1 ( Figure 4B), for which the expression levels were positively correlated with each other (p < 0.01). These three genes may be key genes in the co-morbidity of COVID-19 and ICM.

Construction of PPI Networks to Identify Hub Genes
Individual nodes were hidden to obtain the PPI network. The network was optimized using Cytoscape and 62 nodes were detected ( Figure 4A). The top 10 hub genes in the PPI were screened using four methods (MCC, EPC, Degree, BottleNeck) in the cytoHubba plugin ( Table 1). The intersection of the results identified three pivotal genes, HSP90AA1, HSPA9, and SRSF1 ( Figure 4B), for which the expression levels were positively correlated with each other (p < 0.01). These three genes may be key genes in the co-morbidity of COVID-19 and ICM.

Validation of Hub Genes Using ROC Curve Plots
GSE5406 and GSE164805 were used as internal datasets, and GSE1 GSE211979 were used as external datasets to plot ROC curves for the three HSP90AA1, HSPA9, and SRSF1. All three genes had good diagnostic perform area under the ROC curve (AUC) >70% in both the internal and external datas

Validation of Hub Genes Using ROC Curve Plots
GSE5406 and GSE164805 were used as internal datasets, and GSE116250 and GSE211979 were used as external datasets to plot ROC curves for the three hub genes, HSP90AA1, HSPA9, and SRSF1. All three genes had good diagnostic performances with area under the ROC curve (AUC) >70% in both the internal and external datasets ( Figure 5A-D). We therefore concluded that HSP90AA1, HSPA9, and SRSF1 are pivotal genes in the co-morbidity of COVID-19 and ICM and have considerable diagnostic value. 5A-D). We therefore concluded that HSP90AA1, HSPA9, and SRSF1 are pivotal gene the co-morbidity of COVID-19 and ICM and have considerable diagnostic value.

Construction of TF-mRNA and mRNA-miRNA Regulatory Networks
We constructed a TF-mRNA regulatory network using ENCODE in the NetworkAnalyst tool ( Figure 6). A total of 83 TFs were identified, and 11 of them were associated with only two of the pivotal genes: ZFX, ZNF584, ZNF24, POLR2A, FOXJ2, TAF7, THRAP3, MTA1, GTF2E2, MYC, and ZNF644. The only TF that was associated with all three pivotal genes was KLF9. We also constructed an mRNA-miRNA network using miRTarBase v8.0 in the NetworkAnalyst tool ( Figure 7). A total of 200 miRNAs were identified, and four of them were associated with only two pivotal genes; namely, mir-324-3p, mir-1226-3p, mir-4724-5p, and mir-10a-5p. The only miRNA that was associated with all three hub genes was mir-16-5p. These TFs and miRNAs had high levels of interaction with the three hub genes in the network, indicating they may regulate these genes.

Construction of TF-mRNA and mRNA-miRNA Regulatory Networks
We constructed a TF-mRNA regulatory network using ENCODE in the NetworkAnalyst tool (Figure 6). A total of 83 TFs were identified, and 11 of them were associated with only two of the pivotal genes: ZFX, ZNF584, ZNF24, POLR2A, FOXJ2, TAF7, THRAP3, MTA1, GTF2E2, MYC, and ZNF644. The only TF that was associated with all three pivotal genes was KLF9. We also constructed an mRNA-miRNA network using miRTarBase v8.0 in the NetworkAnalyst tool (Figure 7). A total of 200 miRNAs were identified, and four of them were associated with only two pivotal genes; namely, mir-324-3p, mir-1226-3p, mir-4724-5p, and mir-10a-5p. The only miRNA that was associated with all three hub genes was mir-16-5p. These TFs and miRNAs had high levels of interaction with the three hub genes in the network, indicating they may regulate these genes.

Drug Prediction
We used cMAP to scan the 81 common DEGs for drug-related signatures. A total of 2429 drugs and 171 cMAP classes were predicted. The highest score was for PI3K inhibitor (enrichment score 95.20) and the lowest score was for Tubulin inhibitor (enrichment score −99.48). We screened the top 10 drugs with negative correlations with the common DEGs, and found that 60% of them were Tubulin inhibitors (Table 2), and the highest ranked drugs were vindesine and ON-01910 (both scores −98.77).  Figure 7. miRNA-mRNA regulatory network of the three hub genes. Red indicates hub genes, green indicates miRNAs associated with at least two of the hub genes.

Drug Prediction
We used cMAP to scan the 81 common DEGs for drug-related signatures. A total of 2429 drugs and 171 cMAP classes were predicted. The highest score was for PI3K inhibitor (enrichment score 95.20) and the lowest score was for Tubulin inhibitor (enrichment score −99.48). We screened the top 10 drugs with negative correlations with the common DEGs, and found that 60% of them were Tubulin inhibitors (Table 2), and the highest ranked drugs were vindesine and ON-01910 (both scores −98.77).

Validation by Molecular Docking Simulations
The structure files of vindesine and ON-01910, the drugs with the highest negative scores, were downloaded from PubChem for docking to the three pivotal genes. We used the 3D structure of ON-01910 and the 2D structure of vindesine because its 3D structure was not available in PubChem. The structure files of HSP90AA1, SRSF1, and HSPA9 were downloaded from the RCSB Protein Data Bank (PDB). The results show that vindesine and ON-01910 docked with HSP90AA1 and HSPA9 with binding fractions <−5 for all the complexes (Table 3 and Figure 8). However, for the docking of vindesine and ON-01910 with SRSF1, the ligand conversion failed and docking was not completed. Table 3. Docking of the three hub genes with predicted drug molecules.

Discussion
Cardiovascular complications caused by COVID-19 have been the focus of much attention, including the pathogenesis of COVID-19 comorbidities of cardiovascular diseases, such as acute myocardial infarction and myocarditis [19,20]. Disease signature genes for COVID-19 or ICM have also been reported. However, little research has been focused on ICM and COVID-19 comorbidity [21][22][23][24][25]. We therefore explored potential biomarkers and molecular mechanisms of COVID-19 and ICM.
We first identified 81 common DEGs, 28 up-regulated and 53 down-regulated, between ICM and COVID-19. We hypothesized that these genes were associated with the potential mechanisms of the two co-morbidities. We then performed an enrichment analysis of the DEGs and discovered that the most enriched GO term under biological process was angiogenesis. Amir et al. [26] reviewed pathological angiogenesis caused by endothelial cell dysfunction in COVID-19, and Liao et al. [27] suggested that cell therapy could promote angiogenesis to improve ICM. We therefore speculated that angiogenesis may be one of the phenotypes of the co-morbidities.
To further explore the molecular mechanism of ICM and COVID-19 co-morbidity, we constructed a PPI network of the 81 common DEGs, and found that 62 of them interacted with each other. After screening and validation, we identified three pivotal genes (HSP90AA1, HSPA9, and SRSF1) that may be candidate biomarkers. The internal and external validation suggested that the three pivotal genes may have considerable diagnostic value, which will help further explore COVID-19 and ICM pathogenesis and diagnosis in future studies. These three genes have been shown previously to play crucial roles in angiogenesis.
HSP90AA1 (also known as HSP90A), a therapeutic target for cardiovascular disease and cardiac ageing, has been reported to mediate angiogenesis by direct binding to semisynthetic triterpenoids [28,29], and it was found to be strongly related to the severity degree of COVID-19 disease [30]. HSP90AA1 was found to be a risk factor for the combination of cardiovascular disease and COVID-19 [31]. HSPA9 (also known as Mortalin) can promote angiogenesis and tumor progression by activating the Wnt/β-linked protein signaling pathway [32]. SRSF1 regulates differentiation, proliferation, and adaptive response of immune cells to various stimuli [33]. The YAP-RUNX2-SRSF1 axis has been reported to promote angiogenesis [34], and the circular RNA circSMARCA5 was found to regulate vascular endothelial growth factor A (VEGFA) angiogenesis by binding to SRSF1 [35]. These findings suggest that SRSF1 also has a close association with angiogenesis. Therefore, it is reasonable to speculate that angiogenesis is a key phenotype of ICM and COVID-19.
TFs and miRNAs often act together to regulate gene expression. To identify regulatory factors related to the pivotal genes, we constructed transcription factor and miRNA regulatory networks for these hub genes. One miRNA (hsa-mir-16-5p) and one TF (KLF9), which have been shown to be related to angiogenesis, were found to bind to all three hub genes. For example, knockdown of KLF9 inhibited retinal angiogenesis [36], and knockdown of miR-16-5p increased the level of vascular endothelial growth factor (VEGF) and increased cell viability and angiogenesis [37]. These findings provide insights that can be used to explore the molecular mechanisms underlying the co-pathogenesis of COVID-19 and ICM, and to determine whether hsa-mir-16-5p and KLF9 can regulate angiogenesis by regulating the pivotal genes throughout the course of the diseases.
Finally, drug prediction that targets key genes is feasible. For example, the inhibitor developed for the polymerase PARP-1 and HSPA9 can be used against COVID-19 [38]. We performed drug predictions using cMAP for the 81 common DEGs. To validate their binding effects to the pivotal genes, we selected the drugs with the highest negatively ranked binding scores for molecular docking with the three hub genes. The binding scores of vindesine and ON-01910 were <−5 for docking with HSP90AA1 and HSPA9, which suggested that these two drugs bound spontaneously to the proteins. However, the docking of these two drugs with SRSF1 suggested a ligand conversion error rather than their inability to dock. We considered that the SRSF1 protein structure that we used did not support the conversion of ligands; therefore, whether these drugs can bind to SRSF1 needs to be further investigated. Despite the failure of docking with SRSF1, the docking of vindesine and ON-01910 with HSP90AA1 and HSPA9 tentatively validated the feasibility of these drugs binding to the hub genes, indicating that vindesine and ON-01910 may be disease-specific agents of ICM and COVID-19.
In this study, common DEGs and potential biomarkers of ICM and COVID-19 were discovered for the first time. Our results will help to further clarify the pathophysiology of both these diseases. Our study has some shortcomings. We integrated publicly available sequencing and gene expression data for the analysis, and therefore batch differences in the integration of sequencing chips from different platforms cannot be ruled out. Although external RNA sequencing datasets were used for validation by molecular docking, this was performed only at the mRNA transcriptome level and only with computer simulations. In vivo and ex vivo experimental validation is needed to verify these findings in the future.

Preparation of Data
Gene expression data for many different disorders are publicly available in the GEO (www.ncbi.nlm.nih.gov/geo, accessed on 14 December 2022). GEO2R is a freely accessible internet tool for analyzing differentially expressed genes [39]. The GSE164805 dataset was generated using an Agilent-085982 Arraystar human lncRNA V5 microarray platform. The dataset contains 10 COVID-19 samples and 5 control samples (13 males and 2 females). For all the patients, the diagnosis was confirmed by reverse transcription-polymerase chain reaction (RT-PCR) for COVID-19 levels. Median age was 53.5 and 56 years for healthy individuals and patients with COVID-19, respectively. The GSE5406 dataset was generated using an Affymetrix human genome U133A array platform. The dataset contained 108 ICM samples and 16 control samples. All the patients had New York heart association class 3 to 4 symptoms and left ventricular systolic dysfunction, and all the controls had normal left ventricular function. The external datasets GSE211979 and GSE116250 were both generated by high-throughput sequencing on Illumina HiSeq2500 platforms. GSE211979 contains 16 COVID-19 and 5 control samples. All the patients were diagnosed by RT-PCR for COVID-19. The control samples, which were matched for age and sex, were from the Personalized Environmental and Genetic Study cohort. GSE116250 contains 13 ICM samples and 14 control samples (21 male and 6 female). The patients with ICM were identified based on medical history; all of them had significant obstructive coronary artery disease or myocardial infarction. The control group had no significant cardiac history and had ≥25% shortening of the left ventricular echocardiogram. Median age was 53.5 and 56 years for the control and ICM samples, respectively.

Collection of Common DEGs for COVID-19 and ICM
Differential gene analysis was performed using GEO2R. Statistical significance was defined as adjusted p value < 0.05. The screening thresholds to identify DEGs in the two internal datasets were log 2 (fold change) >0.5 (up-regulated DEGs) and log 2 (fold change) <−0.5 (down-regulated DEGs). The results were collated separately, including the removal of unannotated and/or duplicate genes. The Venn diagrams were integrated to identify DEGs co-regulated or down-regulated by COVID-19 and ICM. The data were visualized using Hiplot (ORG) (https://hiplot.org, accessed on 1 January 2023) [40].

Enrichment Analysis of Common DEGs for COVID-19 and ICM
Functional enrichment analyses were performed on the identified common DEGs using DAVID [41]. Statistical significance was set as p value < 0.05. Bubble plots were used to visualize the GO results, and a string plot was used to visualize the KEGG results.

Construction of PPI Networks to Obtain Hub Genes
Protein functional and regulatory interactions are a major source of complexity in cells. PPI network analysis can help to gain a deeper understanding of disease development. Common DEGs were uploaded to STRING (https://string-db.org, accessed on 5 January 2023) to construct a network of PPIs [42]. The networks were optimized using Cytoscape software 3.8.2 [43]. The top 10 genes in the networks were screened using four methods (MCC, EPC, Degree, BottleNeck) in the cytoHubba plugin to identify hub genes, which were then visualized in Venn diagrams. The expression data of the identified hub genes were integrated. After calculating the mean value for multiple gene probes that corresponded to the expression of the same hub gene, a correlation test was performed on the hub genes.

ROC Curve Validation
GSE5406 and GSE164805 were used as internal datasets, and GSE116250 (ICM) and GSE211979 (COVID-19) were used as external datasets. ROC curves were generated to assess the diagnostic efficacy of the hub genes for ICM and COVID-19.

Construction of TF-mRNA and mRNA-miRNA Regulatory Networks
MiRNAs (miRNAs) are small noncoding RNAs that bind to target mRNAs and regulate their expression by gene silencing or translational repression [44]. TFs also regulate gene expression both spatially and temporally [45]. The ENCODE and miRTar-Base v8.0 databases in the NetworkAnalyst tool [46] (https://www.networkanalyst.ca/ NetworkAnalyst, accessed on 12 January 2023) were used to build regulatory TF-mRNA and mRNA-miRNA networks for the hub genes, respectively. Cytoscape software was used to visualize the two networks.

Prediction of Potential Drugs of Common DEGs
CMap (https://clue.io/, accessed on 2 February 2023) is a database that connects genes, drugs, and disease states through common gene expression profiles [47]. The common COVID-19 and ICM co-expressed DEGs were uploaded to CMap to predict the drugs that may affect these genes. High absolute values of the score are the most relevant; "+" indicates that the trend of the drug effect is the same as the trend of the DEG, and "−" indicates the trend of the drug effect is opposite the trend of the DEG. Drugs with high negative correlation scores are considered the most relevant and can potentially down-regulate up-regulated genes and up-regulate down-regulated genes.

Molecular Docking Validation
Molecular docking is used to explore the binding properties of 3D ligand-receptor structures [48,49] and is widely used in drug discovery. The CB-DOCK2 software for molecular docking uses the tertiary structures of proteins and small ligand molecule. CB-DOCK2 has been used to validate the binding scores of key genes to small drug molecules [50,51]. The 3D structure files of the proteins encoded by the key hub genes were downloaded from RCSB PDB (pdb format) (https://www.rcsb.org/, accessed on 5 February 2023) [52]. The 3D structure files for the predicted drugs were downloaded from PubChem (sdf format) (https://pubchem.ncbi.nlm.nih.gov/, accessed on 5 February 2023). Molecular docking of the hub genes and drugs provide docking scores that indicate how well a drug binds to a protein; low scores indicate good binding.

Conclusions
In this study, HSP90AA1, HSPA9, and SRSF1 were identified as markers for the copathogenesis of COVID-19 and ICM. We also found that the co-pathogenesis of these diseases may be related to angiogenesis. Vindesine and ON-01910 were predicted as potential therapeutic agents for COVID-19 and ICM. These biomarkers and the relationship of COVID-19 and ICM with angiogenesis may help to provide a deeper understanding of the comorbidity of ICM with COVID-19.
Author Contributions: Conceptualization, software, investigation, resources, data curation, writingoriginal draft preparation, S.L., X.Z. and S.T.; methodology, validation, formal analysis, X.X. and W.L.; visualization, Z.Y.; supervision, project administration, funding acquisition, writing-review and editing, W.H. All authors have read and agreed to the published version of the manuscript.
Funding: This research was funded by The National Science Foundation of China, grant number 82074329.
Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The datasets analyzed in this study were acquired from the GEO database (https://www.ncbi.nlm.nih.gov/geo/, accessed on 14 December 2022).