Identification of the Hub Genes and the Signaling Pathways in Human iPSC-Cardiomyocytes Infected by SARS-CoV-2.

Severe Acute Respiratory Syndrome Coronavirus Type 2 (SARS-CoV-2) is an enveloped single-stranded RNA virus that can lead to respiratory symptoms and damage many organs such as heart, kidney, intestine, brain and liver. It has not been clearly documented whether myocardial injury is caused by direct infection of cardiomyocytes, lung injury, or other unknown mechanisms. The gene expression profile of GSE150392 was obtained from the Gene Expression Omnibus (GEO) database. The processing of high-throughput sequencing data and the screening of differentially expressed genes (DEGs) were implemented by R software. The R software was employed to analyze the Gene Ontology (GO) analysis and the Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis. The protein-protein interaction (PPI) network of the DEGs was constructed by the STRING website. The Cytoscape software was applied for the visualization of PPI network and the identification of hub genes. The statistical analysis was performed by the GraphPad Prism software to verify the hub genes. A total of 516 up-regulated genes and 191 down-regulated genes were screened out. The top 1 enrichment items of GO in biological process (BP), Cellular Component (CC), and Molecular Function (MF) were type I interferon signaling pathway, sarcomere, and receptor ligand activity, respectively. The top 10 enrichment pathways, including TNF signaling pathway, were identified by KEGG enrichment analysis. A PPI network was established, consisting of 613 nodes and 3,993 edges. The 12 hub genes were confirmed as statistically significant, which was verified by GSE151879 dataset. In conclusion, the hub genes of human iPSC-cardiomyocytes infected with SARS-CoV-2 were identified through bioinformatics analysis, which may be used as biomarkers for further research.

Introduction

Severe Acute Respiratory Syndrome Coronavirus Type 2 (SARS-CoV-2) is an enveloped single-stranded RNA virus which belongs to the family of coronaviruses and the severe acute respiratory syndrome-related coronavirus (Gorbalenya et al. 2020). As the seventh known coronavirus, it can be transmitted from person to person (Anderson et al. 2020). According to the latest data from Johns Hopkins University in the United States, as of 6:35 a.m. Beijing time, July 23rd, 2020, the cumulative number of confirmed cases and deaths of new coronavirus disease 2019 (COVID-19) worldwide has exceeded 15 million and 620,000, respectively. Currently, the number is still rising rapidly.

The respiratory secretions or droplets of infected people, if inhaled, will cause infection in healthy people. SARS-CoV-2 can spread among people through close contact (Anderson et al. 2020; Shereen et al. 2020; Xu et al. 2020). Fever and cough are frequently seen in patients, but some patients may develop myalgia and fatigue (Xu et al. 2020; Huang et al. 2020). Although the patients mainly show the symptoms of the respiratory system, in fact, the virus can damage many organs at the same time, such as heart, kidney, intestine, brain and liver (Dhakal et al. 2020). It is reported that 23% of severe patients were diagnosed with cardiac injury, a common complication that would exacerbate the disease in its severity and progressing, which includes palpitations or chest pain, myocarditis, acute myocardial infarction (AMI), heart failure (HF), arrhythmia, as well as venous thromboembolic events (VTE) (Yang et al. 2020; Chen et al. 2020; Long et al. 2020; Rizzo et al. 2020). Among them, the AMI and arrhythmia can be seen in patients with COVID-19, and it is more common in ICU patients (Bansal 2020). The receptor-binding domain (RBD) of SARS-CoV-2 might bind tightly with human ACE2, which was strongly expressed in pericytes of human hearts, indicating the potential susceptibility of cardiomyocytes to SARS-CoV-2 (Chen et al. 2020). Additionally, the expression of ACE2 was highly increased in the heart failure patient (Chen et al. 2020). In the direct mechanism, the myocardial injury can be caused by the infiltration of viruses into myocardial tissues, resulting in the death and inflammation of cardiac muscles. In the indirect mechanism, cardiomyocyte stress secondary to pulmonary failure, hypoxemia, and other cardiac inflammation can lead to severe systemic inflammatory response syndrome (SIRS) (Akhmerov and Marbán 2020). The preliminary clinical studies have developed several drugs of therapeutic potential, such as antiviral drugs, nucleoside analogues, neuraminidase inhibitors, therapeutic peptide, RNA synthesis inhibitors, anti-inflammatory drugs, chloroquine phosphate, Lianhuaqingwen capsules and ShuFengJieDu capsules as well (Chakraborty et al. 2020; Lu 2020; Gao et al. 2020).

The high-throughput sequencing technology provides the information of genes, pathophysiological mechanism, diagnosis and therapy, as well as the method of antiviral prevention and vaccination development (Chan et al. 2020). It not only shortens the experimental cycle and reduces the cost of experiments, but also significantly improves the application value of species (such as SARS-CoV-2), the databases of which need to be improved due to the scarcity of their biological information (Pareek et al. 2011).

The exact mechanism of myocardial damage caused by COVID-19 remains unclear. In this study, the gene expression data of human-induced pluripotent stem cell-derived cardiomyocytes (iPSC) infected by SARS-CoV-2 was analyzed through bioinformatics analysis, the result of which is hopefully helpful for studying the mechanism of the infection of myocardium cells by SARS-CoV-2. It is of great significance for the establishment of cardiac antiviral drug screening platform (Bojkova et al. 2020), and helpful to promote the development of clinical targeted drug therapy. The aim of our study is identifying hub genes and signaling pathways that SARS-CoV-2 infected cardiomyocytes, which may be available to screen appropriate and effective biomarkers for further research.

Methods

Data Sources

The gene expression profiles of the GSE150392 dataset was downloaded from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) database of the National Center for Biotechnology Information (NCBI). The database was designed to stockpile primitive sequences, gene expression datasets and platform records. Based on the platform of [GPL18573 Illumina NextSeq 500 (Homo sapiens)], the dataset contains six human iPSC-cardiomyocytes samples, among which the GSM4548303-GSM4548305 served as the infection group while the GSM4548306-GSM4548308 as the control group.

Data Processing

The R software was applied to calibrate and standardize the high-throughput sequencing data of the GSE150392 dataset that was downloaded from the GEO database. The limma package of the R software was employed to screen differentially expressed genes (DEGs) with conditions of |log2FC|≥ 2 and P < 0.05, the pheatmap package to construct the heat-map of DEGs, and the ggplot2 package (Ito and Murphy 2013) to construct the volcano map.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) Analysis of DEGs

The function and pathway enrichment analysis was employed to explain the biological functions of DEGs. The GO analysis, consisting of three categories of cellular component (CC), biological process (BP), molecular function (MF), was applied to describe the characteristics of genes and gene products (Fan and Wei 2020). The KEGG analysis, a public integrated database of biological systems designed for the systematic analysis of gene functions, rendered it possible to understand the pathways which have changed significantly under experimental conditions. The clusterProfiler package in the R software was applied to conduct the GO and KEGG analysis of DEGs, the results of which would be regarded as statistically significant if P < 0.05.

Construction of Protein–Protein Interaction (PPI) Network and Identification and Validation of Hub Genes

The DEGs were uploaded to the STRING website (Szklarczyk et al. 2017) (http://string-db.org/) to establish a PPI network that aimed to analyze the relationship between proteins in DEGs. The visualization and hub-genes screening were performed by Cytoscape software (Su et al. 2014), on which the cytoHubba plug-in was installed to calculate the score of each gene. According to Maximal Clique Centrality (MCC), the 12 genes with the highest degree were identified as hub genes. To confirm the findings from the bioinformatics analysis, the single-cell RNA sequencing’s data of the peripheral blood mononuclear cells (PBMCs) from the GSE150728 dataset and the high-throughput sequencing’s data of human embryonic stem cell-derived cardiomyocytes (hESC-CMs) from the GSE151879 dataset were applied to verify the hub genes. The GraphPad Prism software was used to analyze the data from the GSE150728 dataset and GSE151879 dataset, and P < 0.05 indicated the significant difference.

Construction of gene and miRNA Interactions network

The microRNAs (miRNAs), short non-coding single-stranded RNAs, play a part in pathogenesis of diseases in which they can control a variety of biological functions by regulating gene expression (Chou et al. 2018; Mitash et al. 2020; Saçar Demirci and Adan 2020). After uploading the gene list to the Network Analyst website, the miRTarBase v8.0 database was used to perform network analysis (Chou et al. 2018), which collects information about miRNA-target interactions (MTI) verified by experiment.

Results

Identification of DEGs

Totally 707 DEGs have been screened out according to the criteria of |log2FC|≥ 2 and P < 0.05, including 516 up-regulated genes and 191 down-regulated genes, the results of which were presented in the form of volcano map and heat map (Fig. 1A, B) with the help of ggplot2 package and pheatmap package in the R software.

Fig. 1

Volcano map and heat map of differentially expressed genes (DEGs). A Red dots indicated up-regulated genes and blue dots indicated down-regulated genes. Black dots indicated the rest of the genes with no significant expression change. The threshold was set as followed: P < 0.05 and |log2FC|≥ 2. FC: fold change. B Gene expression data is converted into a data matrix. Each column represents the genetic data of a sample, and each row represents a gene. The color of each cell represents the expression level, and there are references to expression levels in different colors in the upper right corner of the figure

GO Enrichment Analysis of DEGs

The GO enrichment analysis denoted that the parts playing a role of significance in BP section are response to type I interferon, type I interferon signaling pathway, cellular response to type I interferon, etc. (Fig. 2A). As for CC section, sarcomere, contractile fiber part, contractile fiber, I-kappa B/NF-kappa B complex, as well as myofibril rank top 5 (Fig. 2B). Additionally, the MF part presents that the top enrichment items are receptor ligand activity, receptor regulator activity, cytokine receptor-binding, DNA-binding transcription activator activity, RNA polymerase II-specific, and cytokine activity (Fig. 2C). More details of the top 5 significant items in BP, CC, and MF of GO enrichment can be seen in Table 1.

Fig. 2

The GO enrichment analysis and KEGG pathways analysis. BP biological processes, CC cellular components, MF molecular function, KEGG kyoto encyclopedia of genes and genomes; GO gene ontology

Table 1

Significant enrichment of GO terms (top 5 respectively in biological process (BP), cellular component (CC) and molecular function (MF) according to P-value)

Ontology	ID	Description	P-value	Count	Gene ID
BP	GO:0,034,340	Response to type I interferon	2.46627E−22	30	IRF8/SHMT2/IRF1/IRF7/OAS1/MX1/EGR1/OAS3/IFITM1/HLA-E/SP100/ISG20/OAS2/BST2/ISG15/IFIT3/IFNB1/IFI6/SAMHD1/MYD88/CACTIN/IFIT2/USP18/IFIT1/RSAD2/IRAK1/IFI27/MX2/IFI35/STAT2
BP	GO:0,060,337	Type I interferon signaling pathway	9.3015E−22	29	IRF8/IRF1/IRF7/OAS1/MX1/EGR1/OAS3/IFITM1/HLA-E/SP100/ISG20/OAS2/BST2/ISG15/IFIT3/IFNB1/IFI6/SAMHD1/MYD88/CACTIN/IFIT2/USP18/IFIT1/RSAD2/IRAK1/IFI27/MX2/IFI35/STAT2
BP	GO:0,071,357	Cellular response to type I interferon	9.3015E−22	29	IRF8/IRF1/IRF7/OAS1/MX1/EGR1/OAS3/IFITM1/HLA-E/SP100/ISG20/OAS2/BST2/ISG15/IFIT3/IFNB1/IFI6/SAMHD1/MYD88/CACTIN/IFIT2/USP18/IFIT1/RSAD2/IRAK1/IFI27/MX2/IFI35/STAT2
BP	GO:0,009,615	Response to virus	1.03297E−18	47	IL1B/FOSL1/IL6/BIRC3/ZC3H12A/BCL3/PMAIP1/MICB/IRF1/IRF7/OAS1/IL33/MX1/TRIM38/TNFAIP3/OAS3/IFITM1/IL12A/ISG20/OAS2/DTX3L/BST2/F2RL1/IFI44/ISG15/IFIT3/TNFSF4/IFNB1/IFI6/IFI44L/SAMHD1/EIF2AK2/IFIT2/DDX58/DHX58/GBP3/TRIM22/TRIM25/IFIH1/SLFN11/CXCL10/IFIT1/PARP9/RSAD2/IFI27/MX2/STAT2
BP	GO:0,051,607	Defense response to virus	5.53996E−18	39	IL1B/IL6/BIRC3/ZC3H12A/PMAIP1/MICB/IRF1/IRF7/OAS1/IL33/MX1/TRIM38/TNFAIP3/OAS3/ISG20/OAS2/DTX3L/BST2/F2RL1/ISG15/IFNB1/IFI6/IFI44L/SAMHD1/EIF2AK2/DDX58/DHX58/GBP3/TRIM22/TRIM25/IFIH1/SLFN11/CXCL10/IFIT1/PARP9/RSAD2/IFI27/MX2/STAT2
CC	GO:0,030,017	Sarcomere	1.33758E−07	21	ACTC1/SQSTM1/TNNT2/FBXL22/TNNC1/MYL2/MYL4/LMOD2/TCAP/MYL7/SYNPO2L/CSRP3/ITGB1BP2/TNNI3/LDB3/CAV3/ACTA1/TRIM54/HRC/CASQ2/MYZAP
CC	GO:0,044,449	Contractile fiber part	7.06221E−07	21	ACTC1/SQSTM1/TNNT2/FBXL22/TNNC1/MYL2/MYL4/LMOD2/TCAP/MYL7/SYNPO2L/CSRP3/ITGB1BP2/TNNI3/LDB3/CAV3/ACTA1/TRIM54/HRC/CASQ2/MYZAP
CC	GO:0,043,292	Contractile fiber	7.44674E−07	22	ACTC1/SQSTM1/TNNT2/FBXL22/TNNC1/MYL2/MYL4/LMOD2/TCAP/MYL7/SYNPO2L/CSRP3/ITGB1BP2/TNNI3/LDB3/CAV3/ACTA1/TRIM54/CDK5R1/HRC/CASQ2/MYZAP
CC	GO:0,033,256	I-kappaB/NF-kappaB complex	9.59297E−07	4	NFKBIA/BCL3/NFKB2/NFKB1
CC	GO:0,030,016	Myofibril	1.07893E−06	21	ACTC1/SQSTM1/TNNT2/FBXL22/TNNC1/MYL2/MYL4/LMOD2/TCAP/MYL7/SYNPO2L/CSRP3/ITGB1BP2/TNNI3/LDB3/CAV3/ACTA1/TRIM54/HRC/CASQ2/MYZAP
MF	GO:0,048,018	Receptor ligand activity	5.61031E−12	43	IL11/CXCL1/FGF12/GDF15/IL1B/IL6/CCL7/HMGB2/CXCL3/CXCL8/LTA/VEGFC/FGF7/CSF1/INHBA/GREM1/PYY/STC1/CLCF1/TNFSF9/NPPA/LGALS3/OSM/LEFTY2/IL33/CDC42EP2/CXCL6/LIF/CCL26/NRG3/IL12A/TYMP/ADM2/INHBE/TIMP1/GAL/STC2/TNFSF4/IFNB1/BMP6/GHRH/CXCL10/GDF6
MF	GO:0,030,545	Receptor regulator activity	1.44694E−11	44	IL11/CXCL1/FGF12/GDF15/IL1B/IL6/CCL7/HMGB2/CXCL3/CXCL8/LTA/VEGFC/FGF7/CSF1/INHBA/GREM1/PYY/STC1/CLCF1/TNFSF9/NPPA/LGALS3/OSM/LEFTY2/IL33/CDC42EP2/CXCL6/LIF/CCL26/NRG3/IL12A/TYMP/ADM2/INHBE/TIMP1/GAL/STC2/TNFSF4/IFNB1/BMP6/GHRH/CXCL10/GDF6/PXDN
MF	GO:0,005,126	Cytokine receptor-binding	5.59913E−11	29	IL11/IL1B/IL6/CCL7/CXCL3/CCL2/CXCL8/LTA/VEGFC/CSF1/GREM1/CLCF1/TNFSF9/OSM/CXCL6/LIF/CCL26/TRAF1/IL12A/SH2B3/TNFSF4/IFNB1/CCL20/MYD88/SPRED1/ITGA5/CXCL10/RIPK1/PXDN
MF	GO:0,001,228	DNA-binding transcription activator activity, RNA polymerase II-specific	2.3116E−08	36	NR4A2/NR4A1/FOSL1/EGR2/LMO2/HOXA5/MAFF/DBP/DDIT3/JUNB/KLF5/FOSB/FOS/NFATC2/ELF3/SIX4/ETV4/IRF1/NR4A3/MYC/ATF3/NFKB2/IER2/EGR1/CREBRF/CSRNP1/KLF10/ZNF143/JUN/ARNTL/NFKB1/ETS1/BACH1/EPAS1/STAG1/ZEB2
MF	GO:0,005,125	Cytokine activity	3.21362E−08	21	IL1B/IL6/CCL7/CXCL3/CXCL8/LTA/CSF1/INHBA/GREM1/CLCF1/TNFSF9/OSM/IL33/CXCL6/LIF/CCL26/IL12A/TIMP1/TNFSF4/IFNB1/CXCL10

KEGG Pathway Enrichment Analysis of DEGs

The biological functions of DEGs were identified through the KEGG pathway enrichment analysis. The result of the analysis of is shown in Fig. 2D. According to the P-value, the top 10 enrichment pathways were obtained, including TNF signaling pathway, NF–kappa B signaling pathway, Legionellosis, Cytokine–cytokine receptor interaction, IL-17 signaling pathway, Measles, Influenza A, Epstein-Barr virus infection, NOD-like receptor signaling pathway and Pertussis, as shown in Table 2.

Table 2

The top 10 significantly enriched KEGG pathways (according to P-value)

ID	Description	P-value	Count	Gene ID
hsa04668	TNF signaling pathway	1.63058E−17	30	CXCL2/NFKBIA/CXCL1/IL1B/IL6/BIRC3/BCL3/PTGS2/CXCL3/CCL2/JUNB/LTA/FOS/VEGFC/CSF1/MAP3K8/IRF1/SOCS3/CXCL6/TNFAIP3/LIF/TRAF1/ICAM1/MLKL/IFNB1/CCL20/JUN/CXCL10/NFKB1/RIPK1
hsa04064	NF-kappa B signaling pathway	7.05806E−13	24	CXCL2/NFKBIA/CXCL1/IL1B/BIRC3/PTGS2/CXCL3/RELB/GADD45B/CXCL8/LTA/NFKB2/TNFAIP3/TRAF1/ICAM1/TLR4/MYD88/DDX58/TRIM25/IL1R1/CYLD/NFKB1/IRAK1/RIPK1
hsa05134	Legionellosis	2.33394E−10	16	CXCL2/NFKBIA/CXCL1/IL1B/IL6/CXCL3/CXCL8/NFKB2/CLK1/IL12A/HSPA6/TLR4/CLK4/MYD88/NFKB1/C3
hsa04060	Cytokine-cytokine receptor interaction	2.36757E−10	37	IL11/CXCL2/CXCL1/GDF15/IL1B/IL6/CCL7/CXCL3/CCL2/CXCL8/LTA/CSF1/INHBA/CLCF1/TNFSF9/OSM/IL33/CXCL6/LIF/CCL26/IL15RA/IL12A/INHBE/TNFRSF6B/IL1RL1/TNFSF4/IFNB1/CCL20/BMP6/OSMR/IL1R1/CXCL10/GDF6/IL13RA1/TNFRSF10B/IL4R/TNFRSF10D
hsa04657	IL-17 signaling pathway	3.52987E−10	20	CXCL2/NFKBIA/CXCL1/IL1B/FOSL1/IL6/CCL7/PTGS2/CXCL3/CCL2/CXCL8/FOSB/FOS/CXCL6/TNFAIP3/S100A8/CCL20/JUN/CXCL10/NFKB1
hsa05162	Measles	6.21493E−10	24	NFKBIA/IL1B/IL6/FOS/IRF7/OAS1/MX1/TNFAIP3/OAS3/RAB9B/IL12A/HSPA6/TLR4/OAS2/IFNB1/JUN/EIF2AK2/MYD88/DDX58/IFIH1/NFKB1/IRAK1/EIF2AK3/STAT2
hsa05164	Influenza A	2.11084E−09	26	NFKBIA/IL1B/IL6/CCL2/CXCL8/IRF7/OAS1/SOCS3/IL33/MX1/OAS3/ICAM1/IL12A/TLR4/OAS2/IFNB1/EIF2AK2/MYD88/DDX58/TRIM25/IFIH1/CXCL10/TNFRSF10B/RSAD2/NFKB1/STAT2
hsa05169	Epstein-Barr virus infection	1.84108E−08	27	NFKBIA/IL6/HES1/RELB/GADD45B/MYC/NFKB2/IRF7/OAS1/TNFAIP3/OAS3/ICAM1/HLA-E/OAS2/ISG15/CD44/IFNB1/JUN/EIF2AK2/MYD88/DDX58/TAP1/CXCL10/NFKB1/IRAK1/RIPK1/STAT2
hsa04621	NOD-like receptor signaling pathway	2.57364E−08	25	CXCL2/NFKBIA/CXCL1/IL1B/IL6/BIRC3/CXCL3/CCL2/CXCL8/CASP4/IRF7/OAS1/TNFAIP3/OAS3/TLR4/OAS2/PLCB2/IFNB1/JUN/MYD88/GBP3/PANX1/NFKB1/RIPK1/STAT2
hsa05133	Pertussis	1.99721E−07	15	IL1B/IL6/IRF8/CXCL8/FOS/IRF1/CXCL6/IL12A/TLR4/JUN/MYD88/ITGA5/NFKB1/IRAK1/C3

Construction of PPI Network and Identification and Validation of Hub Genes

A PPI network was constructed in the STRING website, which consists of 82 nodes and 213 edges. By visualizing and screening hub genes with Cytoscape software, a PPI network and top 12 genes were obtained and presented in Fig. 3A, B. The top 12 genes were: ISG15, MX1, IFIT1, IFIT2, RSAD2, OAS2, IFIT3, MX2, OAS1, OAS3, IFI6 and DDX58, which were confirmed as statistically significant verified by the GSE151879 dataset (Fig. 4A). The relative expression levels of 10 hub genes, namely ISG15, MX1, IFIT1, RSAD2, OAS2, IFIT3, MX2, OAS1, OAS3 and IFI6, were also consistent with single-cell RNA sequencing of PBMCs from the GSE150728 dataset, while DDX58 and IFIT2 showed no statistically significant difference (Fig. 4B).

Fig. 3

The PPI network of DEGs and the top 12 hub genes of PPI network

Fig. 4

A Verification of hub genes with GSE151879 dataset. P-value < 0.05 is considered to be statistically significant. B Verification of hub genes with GSE150728 dataset. P-value < 0.05 is considered to be statistically significant

Construction of gene and miRNA Interactions network

After network visualization, Gene-miRNA interactions were presented in Fig. 5. According to the criterion of the degree, the miRNA-hsa-mir-335-5p was finally screened out.

Fig. 5

The gene-miRNA network of DEGs based on miRTarBase v8.0 database

Discussion

In this study, the gene expression profile of the GSE150392 dataset was downloaded from the GEO database, and 707 DEGs were identified, including 516 up-regulated DEGs and 191 down-regulated DEGs. In addition, the GO enrichment analysis and KEGG pathway analysis were applied to analyze the functional enrichment of DEGs. A PPI network was applied to identify the hub genes that might play a significant role in the infection of cardiomyocytes by SARS-CoV-2.

As presented in the figure that described the results of GO enrichment analysis of DEGs, it is believed that the smaller the P-value is, the more significant the extent of enrichment in the GO analysis is. As for biological processes (BP), DEGs are most significantly enriched in response to type I interferon, type I interferon signaling pathway and cellular response to type I interferon. Type I interferon has a broad-spectrum antiviral feature against RNA viruses, which may lead to innate and adaptive immune response to function against viruses. Besides, IFN-α promotes the expression of ACE2 that is strongly expressed in human hearts. The RBD of SARS-CoV-2 binding ACE2 leads to infection, while SARS-CoV-2 can employ species-specific interferon to drive up-regulation of ACE2 to promote infection (Chen et al. 2020; Ziegler et al. 2020; Mantlo et al. 2020). In terms of cellular components (CC), sarcomere and contractile fiber play a role of significance in the GO enrichment. The direct damage of cardiomyocytes, cardiac inflammation, immune response, myocardial interstitial fibrosis, and hypoxia caused by SARS-CoV-2 might act on sarcomere and contractile fiber probably, all of which may result in functional changes in the sarcomere and contractile fiber, and gradually evolve into myocardial injury (Babapoor-Farrokhran et al. 2020). In the respect of molecular function (MF), receptor ligand activity and receptor regulator activity are most significantly enriched in DEGs. The ACE2 appears to be the crucial functional receptor of SARS-CoV-2 (Coutard et al. 2020; Walls et al. 2020). Interaction between viral protein ligand and cell receptor of SARS-CoV-2, a critical step in replication of viruses, is also a crucial process of virus infection, which might be determined by the type of chemical interactions, activity between receptor and ligand, and the regulator activity (Ortega et al. 2020).

The KEGG pathway enrichment analysis discovered that quite a few significantly enriched KEGG pathways rated by P-value, such as TNF signaling pathway, NF-kappa B (NF-κB) signaling pathway, and Legionellosis, are highly relevant to the infection of SARS-CoV-2. As for the TNF signaling pathway, TNF is the crucial inflammatory cytokine resulting in various acute and chronic inflammations (Soy et al. 2020). TNF-α level was observed much higher in serum of patients with COVID-19, implying that it is positively correlated with the severity of disease (Soy et al. 2020). With the progress of SARS-CoV-2 infection, monocytes and macrophages increased, resulting in cytokine storm, which meant the release of cytokine and pro-inflammatory cytokines, including IL-1, IL-8, and TNF (Soy et al. 2020; Runfeng et al. 2020). The inflammatory response triggered by the SARS-CoV-2 in myocardium drives an up-regulation of TNF signaling pathway. In terms of the following pathway of NF-kappa B (NF-κB), it is involved in inflammation, innate and adaptive immune response, pathological development of tumor or cancer, and even the critical regulators of stress responses, apoptosis and differentiation (Baldwin 2001; DiDonato et al. 2012; Oeckinghaus et al. 2011). SARS-CoV-2 triggered the stimulation of humoral and cellular immune response as well as mitogen-activated protein kinase (MAPK) pathway and NF-kappa-B signaling pathway, and gradually caused deregulation of immune system and up-regulation of inflammatory pathways, especially NF-kappa B (Mozafari et al. 2020; Wu and Yang 2020; Li et al. 2020). The higher expression of NF-kappa B conducted the induction of the naïve T cell activation and proliferation, promoted the expression of chemokine expression and other immune response, and finally resulted in inflammatory cytokine storm in myocardium (Mozafari et al. 2020; Hoesel and Schmid 2013).

The 12 hub genes of cardiomyocytes specifically infected by SARS-COV-2 were identified, namely ISG15, MX1, IFIT1, IFIT2, RSAD2, OAS2, IFIT3, MX2, OAS1, OAS3, IFI6 and DDX58, which may become the biomarkers. Among the 12 genes, Anna et al. (Rahnefeld et al. 2014) found that ISG15, as a part of innate immunity of cardiomyocytes, plays a significant role in the inhibition of viral replication in mouse cardiomyocytes infected by Coxsackievirus B3 (CVB3). Besides, SARS-CoV-2 can induce a large number of ISG expressions, either as effectors, regulators or both, which is the antiviral response of the body (Clemente et al. 2020). It is suggested that ISG15 may play a crucial role in cardiomyocytes infected by SARS-CoV-2. As for another hub gene (MX1), Eshwar et al. have shown that the β4 subunit of Cav1.2 channels can promote the expression of IFN-related genes (including MX1) in cardiomyocytes, thus alleviating viral infection (Tammineni et al. 2018). IFIT1, IFIT2, IFIT3 and IFIT5 constitute the human IFIT gene family, which expressed very low in most cell types but can be significantly increased by interferon therapy, viral infection and pathogen-related molecular pattern (PAMP) (Pidugu et al. 2019). SARS-CoV-2 can induce strong antiviral response by up-regulating antiviral factors such as OAS1-3, IFIT1-3 and CXCL9 / 10 / 11 of Th1 chemokine, and reduce ribosomal protein transcription (Lieberman et al. 2020).

Hsa-miR-335-5p regulates significant transcription factors and participates in activating the pathways affiliated with inflammation and angiogenesis. It may activate WNT and TGF-β signaling pathways, both of which exert pleiotropic and multifunctional effects, regulating a considerable number of biological processes, TGF-β and WNT signaling pathways (Kay et al. 2019; Esquinas et al. 2017). In the pathogenesis of cardiac remodeling and fibrosis, TGF-β plays an important role (Dobaczewski et al. 2011).

The gene expression profiles of normal human iPSC cardiomyocyte samples and human iPSC cardiomyocyte samples infected by SARS-CoV-2 were studied through bioinformatics analysis, and finally, the hub genes of human iPSC-cardiomyocytes after SARS-CoV-2 infection were identified. Although sequencing data from other datasets were used for verification, further experimental research analysis (such as qRT-PCR) still need to be conducted to acquire abundant data to verify the predicted results obtained from the bioinformatics analysis. There are still some limitations to this study. It failed to understand clearly the mechanism of several dominant genes being involved in the infection of cardiomyocytes by SARS-CoV-2. The other limitation is that the small sample size and data sets of this study would affect the accuracy of the conclusions. Further research is entailed to demonstrate the molecular mechanism of these biological events involved in the infection of cardiomyocytes by the novel coronavirus.

Conclusion

In conclusion, the RNA-Seq data of human iPSC-cardiomyocytes infected by SARS-CoV-2 was studied through bioinformatics analysis. Finally, the hub genes associated with SARS-CoV-2 infection were identified, which may reveal the mechanism of heart damage caused by SARS-CoV-2 or suggest biomarkers.