Inflammatory activation and immune cell infiltration are main biological characteristics of SARS-CoV-2 infected myocardium

ABSTRACT The severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2) can target cardiomyocytes (CMs) to directly invade the heart resulting in high mortality. This study aims to explore the biological characteristics of SARS-CoV-2 infected myocardium based on omics by collecting transcriptome data and analyzing them with a series of bioinformatics tools. Totally, 86 differentially expressed genes (DEGs) were discovered in SARS-CoV-2 infected CMs, and 15 miRNAs were discovered to target 60 genes. Functional enrichment analysis indicated that these DEGs were mainly enriched in the inflammatory signaling pathway. After the protein-protein interaction (PPI) network was constructed, several genes including CCL2 and CXCL8 were regarded as the hub genes. SRC inhibitor saracatinib was predicted to potentially act against the cardiac dysfunction induced by SARS-CoV-2. Among the 86 DEGs, 28 were validated to be dysregulated in SARS-CoV-2 infected hearts. Gene Set Enrichment Analysis (GSEA) analysis of Kyoto Encyclopedia of Genes and Genomes (KEGG) showed that malaria, IL-17 signaling pathway, and complement and coagulation cascades were significantly enriched. Immune infiltration analysis indicated that ‘naive B cells’ was significantly increased in the SARS-CoV-2 infected heart. The above results may help to improve the prognosis of patients with COVID-19.


Introduction
Coronavirus disease 2019 (COVID- 19), which results from the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), is spreading all over the world and thus causing a serious global health crisis [1]. SARS-CoV-2 has a high affinity to angiotensin-converting enzyme 2 (ACE2) and ACE2 plays a pivotal role in the pathogenesis of this infectious disease. ACE2 is ubiquitous and widely expressed in human tissues especially in alveolar cells [2]. As expected, the lung is one of the most important target organs during the SARS-CoV-2 attack on the body. It has been reported that some patients progress to severe respiratory failure and acute respiratory distress syndrome (ARDS) which requires mechanical ventilation, and that the mortality rate of these patients is as high as 20-40% [3].
Except for the lung, ACE2 is also highly expressed in the cardiovascular system [4]. In fact, in the systemic symptoms caused by SARS-CoV-2 infection, myocardial injury is one of the important pathogenic features which further increases the risk of mortality [5]. According to previous reports on SARS-CoV, the existence of viral RNA and macrophage infiltration could be simultaneously detected in clinical heart samples, suggesting direct and indirect viral-induced damage to the myocardium [6]. Perez-Bermejo JA et al [7] performed experiments and discovered that SARS-CoV-2 was more inclined to infect cardiomyocytes (CMs) over fibroblasts and endothelial cells. However, the downstream effects on cardiomyocytes after infection with this virus are not delineated.
Therefore, in this study, the differentially expressed genes (DEGs) in SARS-CoV-2 infected three types of human CMs were re-analyzed, and the overlapping genes were enrolled for analysis of protein-protein interaction (PPI), Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), as well as their upstream regulators. Besides, the potential therapeutic drugs are predicted. Furthermore, these DEGs in SARS-CoV-2 infected human hearts were re-analyzed and used to validate the DEGs of in-vitro results. Finally, based on the in vivo results, Gene Set Enrichment Analysis (GSEA) analysis was performed and immune infiltration was evaluated. We hope that the above bioinformatics analysis will show the biological characteristics during SARS-CoV-2 infected myocardium.

Bioinformatics datasets
Three datasets GSE150392, GSE151879 and GSE169241 were obtained from Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) database. In GSE150392, human induced pluripotent stem cell-derived cardiomyocytes (hiPSC-CMs) were infected with Mock or SARS-CoV-2. In GSE151879, human embryonic stem cells-derived cardiomyocytes (hESC-CMs) and adult human cardiomyocytes (AHC Ms) were both infected with Mock or SARS-CoV-2. In GSE169241, heart samples of three COVID-19 patients and five controls were collected. All the gene expression profiles of the aforementioned datasets were detected by high-throughput sequencing.

Identification of upstream TFs and miRNAs
The expression of DEGs could be specifically regulated by TFs at transcriptional level [13] and miRNAs at post-transcriptionally [14]. The transcription factors (TFs) of these DEGs were predicted with FunRich (http://www.funrich.org/) software [15], and the top 10 TFs were picked out. The upstream miRNAs were searched on miRNET website (https://www.mirnet.ca/) [16], and the miRNA-mRNA network was subsequently constructed using the website.

Functional enrichment analysis
GSEA is a knowledge-based approach for analyzing groups of genes that share common biological functions, chromosomal location, or regulation [17]. GO analysis mainly reflects three independent ontologies of group genes: molecular function (MF), biological process (BP) and cellular component (CC) [18]. WEB-based Gene SeT AnaLysis Toolkit (WebGestalt, http://www.webgestalt.org/), a userfriendly functional enrichment analysis webtool [19], was applied for GSEA and GO analysis. KEGG database provided biological pathway information of DEGs enrichment [20]. Enrichr website (https:// amp.pharm.mssm.edu/Enrichr/), a comprehensive, freely available gene set enrichment analysis web server [21], was used to explore pathway information.

PPI network construction and analysis
The 'Search Tool for Retrieval of Interacting Genes/ Proteins' database (STRING, https://string-db.org/) integrates all available information regarding protein-protein interaction [22], and all DEGs were mapped into the database to construct the complete PPI network. Subsequently, the network was imported into Cytoscape software (https://cytos cape.org/), an open-source software platform for visualizing complex networks [23], to screen out hub nodes with plug-in 'cytoHubba' according to the indexes Degree, Betweenness, and Closeness. Furthermore, the plug-in 'MCODE' was used to extract the clustering function modules. The criteria used was as follows: degree cutoff = 2, node density cutoff = 0.1, node score cutoff = 0.2, k-core = 2 and max depth = 100.

Potential therapeutic drugs prediction
The Connectivity Map (cMap) source can provide information of connection between gene expression signatures and small molecules, genes or disease [24]. In order to search for drugs against virus-induced myocardial injury, drug repurposing analysis was performed on the cMap database by using the online application CLUE (https://clue.io) [25]. The first 10 molecules with the lowest scores were collected.

Immune infiltration analysis
CIBERSORT (https://cibersort.stanford.edu/) is an analytical online webtool that provides an estimation of the abundance of member cell types in a mixed cell population by using gene expression data [26]. The immune infiltration of 22 kinds of immunocytes was evaluated with CIBERSORT and analyzed with the wilcoxon test method to compare the difference between infected and uninfected hearts. A p-value < 0.05 was considered statistically significant. Pearson correlation between infiltrating immunocytes was calculated in GraphPad Prism 8.0.

Results
In the present research, we aimed to explore the biological characteristics of SARS-CoV-2 infected myocardium by bioinformatic tools. After reanalyzing the datasets of SARS-CoV-2 infected CMs, the DEGs were picked out and enrolled for GO, KEGG analysis, and PPI network construction. The upstream regulator including TFs and miRNAs that are potentially responsible for the expression of DEGs were screened. In order to find possible therapeutic drugs, the DEGs were collected for analysis in the specific database. Finally, the dataset of the hearts of COVID-19 patients was re-analyzed to confirm the in-vitro results, and further used for detecting the immune infiltration.
In order to more accurately screen out the DEGs after the virus interfered with CMs, we took the intersection of up-and down-regulated DEGs among three groups. The results showed that there were 62 genes up-regulated in all three groups ( Figure 1(d), Table 1), and only 24 genes were down-regulated in all three groups (Figure 1(e), Table 1). In short, these 86 DEGs were considered to be the genes that are most likely differentially expressed after SARS-CoV-2 infected the CMs.
Besides, the upstream miRNAs were predicted and 15 miRNAs were found to potentially regulate the expression of 60 of 86 DEGs, and they formed a complex network containing 132 edges ( Figure 2 (b), Supplementary Table 1). Among all the miRNAs, miR-1-3p and miR-27a-3p had the most target genes.

Analyze GO and KEGG
We further assessed the biological role of 86 DEGs with GO analysis. The results showed that these DEGs were enriched in response to stimulus, biological regulation, and metabolic process under BP categories (Figure 3(a)). In terms of MF, the functions of DEGs were mainly involved in protein binding, ion binding and nucleic acid binding (Figure 3(a)). These DEGs were distributed in different structures of cells including membrane, nucleus, and extracellular space (Figure 3(a)).
Biological pathway was analyzed with KEGG database [20]. As shown (Figure 3(b)), several signaling pathways were significantly enriched including TNF signaling pathway, Cytokinecytokine receptor interaction, IL-17 signaling pathway, NF-kappa B signaling pathway, etc., indicating that the function of DEGs was closely related to inflammatory phenotype.

Construction of the PPI network and hub genes analysis
To construct the PPI network, we mapped all the 86 DEGs into STRING database and generated a complex network structure, in which 54 nodes and 121 edges were retained (Figure 4(a)). After 'MCODE' [27] analysis, we revealed two modules where the first one scored 8.75 (Figure 4(b)). 'cytoHubba' [28] was another way to evaluate the nodes through 'degree', 'betweenness' and 'closeness'. The top 10 genes under three dimensions were displayed (Figure 4(c-e)). Altogether, several genes such as CCL2, CXCL3, CXCL8, CXCL2, IRF1, etc. were potential hub ones.

Validation of the DEGs expression in the dataset of SARS-CoV-2 infected heart
To fully understand the biological characteristics of virus-infected myocardium, we re-analyzed the DEGs between COVID-19 infected hearts and uninfected hearts. As indicated ( Figure 5(a)), there were 1765 upregulated and 1509 downregulated DEGs in the SARS-CoV-2 infected heart. More than that, 15 upregulated and 13 downregulated genes of the 86 DEGs were also dysregulated in heart samples ( Figure 5(b)), showing that the in vitro SARS-CoV-2 infected CMs could simulate in vivo infection in the heart to some extent. All the 28 DEGs were ordered according to the size of their logFC in GSE169241 ( Figure 5(c)), and CCL2 was found to be the most increased while OXCT1 decreased the most in SARS-CoV-2 infected hearts. We also performed GSEA analysis to uncover the biological pathways involved. As shown, 31 KEGG terms including 'Malaria', 'IL-17 signaling pathway', 'Complement and coagulation cascades', etc. were  positively associated with the pathological state, while 25 terms including 'propanoate metabolism', 'valine, leucine, and isoleucine degradation', and 'Citrate cycle (TCA cycle)', etc. were negatively associated with the disease progression ( Figure 6  (a,b)).

Assessing the immune infiltration in the SARS-CoV-2 infected heart
The infection of SARS-CoV-2 was able to induce abnormal response of the inflammatory signal pathway. We, therefore, assessed the immune infiltration of SARS-CoV-2 infected hearts (Figure 7 (a)). The results showed that there were closely positive correlations and negative correlations between the infiltration fraction of multiple inflammatory cells (Figure 7(b)). Furthermore, the proportion of 'naïve B cells' was significantly upregulated (P = 0.036) after the SARS-CoV-2 infection (Figure 7(c)), while 'Macrophages M2 and 'Mast cells resting' showed a downward trend (P = 0.071 for both). This provides a clue that naive B cells may be the focus of intervention.

Discussion
There is no doubt that SARS-CoV-2 induced COVID-19 mainly focused on the respiratory system and results in severe respiratory failure ultimately causing death [29]. Besides, one more important reason for the death of COVID-19 patients is the destruction of the cardiovascular system. Myocardial injury is common during SARS-CoV-2 infection, and researchers have found that it is significantly associated with a higher risk of mortality [30,31]. Therefore, it is urgent to explore the biological effects of SARS-CoV-2 on the heart.  Genes are key participants in life activities and disease evolution. In this study, we firstly collected three independent datasets from the GEO database and re-analyzed them. These three datasets used SARS-CoV-2 or mock treatment to interrogate a response from cardiomyocytes originating from hiPSC, HESC, and adult humans, and detected via high throughput sequencing. The final results  revealed thousands of differentially expressed genes in each dataset. However, only 86 DEGs were concurrently up-regulated or downregulated in all datasets, indicating the existence of several internal biological differences among the above CMs from the three sources when exposed to SARS-CoV-2. Therefore, these overlapped genes were selected as the most potential abnormally expressed factors in SARS-CoV-2 infected CMs. To verify the reliability of these genes, we reanalyzed the dataset of heart samples from clinical COVID-19 patients and uninfected individuals and found that 28 genes of the 86 DEGs were also significantly differentially expressed in vivo, suggesting that the in vitro model was reliable to some extent.
To construct the integrated regulatory network, we predicted the upstream transcription factors and miRNAs with bioinformatic tools. In terms of the transcription factors listed, only NFIC was significantly enriched while others showed no significance. MiRNAs are tissue-specific regulators that can inhibit gene expression. We further predicted the possible miRNAs targeting the 86 DEGs in muscle to construct a precise network. We found 15 miRNAs that potentially regulate a total of 60 genes from 86 DEGs and miR-1-3p is most widely associated with these DEGs. According to previous research, miR-1-3p was increased in peripheral blood and positively correlated with the level of myocardial damage in acute viral myocarditis [32] and therefore whether its circulating form could contribute to diagnose and even judge the prognosis in SARS-CoV-2 induced myocardial injury is worth exploring. Besides, its extensive association with multiple DEGs suggests that it is also a potentially valuable intervention target.
There is still insufficient understanding of this virus-induced myocardial disorder. Based on all dysregulated genes, we evaluated the characteristics of the above biological dysfunction through GO and KEGG analysis. In terms of the biological process in GO, we found that response to stimulus, biological regulation, and metabolic process were significantly enriched. The analysis of biological pathways by KEGG provides abundant information. As indicated in the analysis, these DEGs were enriched in 'TNF signaling pathway', 'Cytokinecytokine receptor interaction', 'IL-17 signaling pathway', and 'NF-kappa B signaling pathway', whose roles were also vital in other viral myocarditis [33,34]. This suggested that antiinflammatory strategies will likely help to combat myocardial injury in COVID-19 patients. We also focused on the KEGG signaling pathway by GSEA analysis with the results of in vivo dataset, and confirmed the importance of 'IL-17 signaling pathway' in the disease progression once again. In fact, activation of IL-17 signaling pathway, to a certain extent, is a landmark event of SARS-CoV-2 infection [35], and mediated the cytokine storm caused by SARS-CoV-2 infection [36]. More than that, 'Complement and coagulation cascades' were proven to actively participate in this cardiac anomaly, which highly suggests that myocardial injury in COVID-19 may be partly from microcirculatory thrombosis.
To identify the key genes responsible for exacerbating SRARS-CoV-2 damage to CMs, we constructed a PPI network with all 86 DEGs and further analyzed the nodes by using 'MCODE' and 'cytohubba' in cytoscape software. Several genes including CCL2, CXCL8, PTGS2, and IRF1 were considered to be hub genes because they were screened out according to different standards simultaneously. Chemokine (C-C-motif) ligand 2 (CCL2), also known as monocyte chemoattractant protein-1 (MCP-1), is a member of the chemokine family, and responsible for attracting leukocytes to sites of infection or injury to mediate defense and repair [37]. In a recent research, SARS-CoV-2 infected CMs could induce CCL2 secretion and further recruit monocyte [9][], indicated a vital role of CCL2 in the pathological mechanism. CXCL8, also called IL-8, is a prototypical member of the CXC family, and is usually responsible for the recruitment and activation of inflammatory cells through interaction with cell surface receptor GPCR, CXCR1 and CXCR2 [38]. In COVID-19, CXCL8 expression was higher only in severe patients, and positively correlated with the percentage of neutrophils [39], which proposed that CXCL8 leads to deterioration of the condition by recruiting neutrophils. CXCL8, along with IL-37 and CRP, could be combined into a highly sensitive model to effectively differentiate severe cases from moderate ones in the COVID-19 population [40]. Therefore, the above hub genes such as CCL2 and CXCL8 may be key participants and diagnostic indicators in COVID-19 patients who experienced myocardial injury.
Due to the fact that these hub genes were closely associated with inflammatory cells in function, we also evaluated the immune infiltration in the injured heart sample with the in vivo dataset. After the execution in the CIBERSORT webtool and subsequent statistical analysis, only the proportion of 'naive B cells' was significantly increased in SARS-CoV-2 infected hearts, while 'Macrophages M2ʹ and 'Mast cells resting' decreased with parameter p-value close to 0.05. The naive B lymphocytes are a major mediator of adaptive immunity in mammals [41]. In a recent research by Wu et al [42]. the immune infiltration landscape in the lung tissues of COVID-19 patients was estimated, and naive B cells infiltration was found to be high and may be the main cause of the over-active humoral immune response. Furthermore, it was found that the increased B cell infiltration was associated with poorer clinical outcomes, and may act as a trigger of severe respiratory and pulmonary symptoms of COVID-19. Several specific molecules that can target B cells have been shown to be beneficial against the COVID-19 pandemic [43,44]. Therefore, the increased B cell infiltration is a potential culprit of functional deterioration in the hearts of COVID-19, and appropriate inhibition of naive B cells might be an effective strategy to alleviate COVID-19-induced heart dysfunction/ damage.
Exploring therapeutic drugs is of great significance to improve the prognosis of COVID-19 patients. We matched 86 DEGs into the cMAP database and obtained scores of all 2429 small molecules. A series of compounds were identified to potentially reduce the SARS-CoV-2 induced CMs injury. The antiviral saracatinib, one of the most likely therapeutic drugs, has been considered to be effective against SARS-CoV-2 in two different research models [45,46]. However, whether saracatinib or other predicted small molecules could prevent myocardial tissue damage caused by SARS-CoV-2 remains to be confirmed.

Conclusion
In conclusion, our study identified potential biomarkers in SARS-CoV-2 infected myocardium, and found potential immune cell intervention targets and valuable therapeutic drugs. However, further studies are needed to confirm our preliminary results.

Highlights
• The miRNA-mRNA network and PPI network will improve our understanding of the disease development and can be applied for exploring potential disease biomarkers. • Fraction of 'naive B cells' was significantly increased in SARS-CoV-2 infected myocardium, and may be a potential intervention target. • Several small moleculars especially saracatinib were believed to fight disease progression.