Targeting hub genes and pathways of innate immune response in COVID-19: A network biology perspective.

The current pandemic of 2019 novel coronavirus disease (COVID-19) caused by a novel virus strain, 2019-nCoV/SARS-CoV-2 have posed a serious threat to global public health and economy. It is largely unknown how the human immune system responds to this infection. A better understanding of the immune response to SARS-CoV-2 will be important to develop therapeutics against COVID-19. Here, we have used transcriptomic profile of human alveolar adenocarcinoma cells (A549) infected with SARS-CoV-2 and employed a network biology approach to generate human-virus interactome. Network topological analysis discovers 15 SARS-CoV-2 targets, which belongs to a subset of interferon (IFN) stimulated genes (ISGs). These ISGs (IFIT1, IFITM1, IRF7, ISG15, MX1, and OAS2) can be considered as potential candidates for drug targets in the treatments of COVID-19. We have identified significant interaction between ISGs and TLR3 agonists, like poly I: C, and imiquimod, and suggests that TLR3 agonists can be considered as a potential drug for drug repurposing in COVID-19. Our network centric analysis suggests that moderating the innate immune response is a valuable approach to target COVID-19.

Introduction

The current pandemic of Coronavirus Disease-2019 (COVID-19) has led to over 5 million confirmed cases and almost 3.5 lakhs fatalities worldwide, since its emergence in late 2019. COVID-19 is caused by a novel virus coronavirus strain, SARS-CoV-2, an enveloped, positive-sense, single-stranded RNA, β-coronavirus of the family Coronaviridae [1].

Nearly 80% of COVID-19 patients show mild symptoms, such as cough and fever, and do not require hospitalization, but the remaining 20% do require [2]. However, of those 20%, almost half of them develop severe respiratory failure in the form of fatal acute respiratory distress syndrome (ARDS) [3]. Along with this fatal form, dysregulated immune response also mediated the severe COVID-19 pathogenesis. The immune dysregulation, called hypercytokinemia or “cytokine storm,” is frequently associated with ARDS [4]. Till date, it remains unclear how SARS-CoV-2 disrupt the host innate immune response. Several recent studies have reported the dysregulated secretion of proinflammatory cytokines in COVID-19 [[5], [6], [7]]. Recently, Melo et al. [5] reported a moderate interferon (IFN) response to SARS-CoV-2 infection in primary cells and showed that IFN can reduce SARS-CoV-2 replication in vitro. Zhou et al. [7] and Xiong et al. [6] recently examined the innate immune response to SARS-CoV-2 infection in bronchoalveolar lavage fluid (BALF). The differentially expressed genes were enriched in inflammatory pathways including Chemokine Signaling [7]. Whereas, in studies by Xiong et al. [6], the upregulated genes were largely related to viral infection with the most enriched biological processes being protein targeting to membrane and ER. These studies reported the significant upregulation of a subset of interferon-stimulated genes (ISGs) which are directly related to antiviral activity, such as ISG15, IFIH1, MX1, OAS1-3, IFITMs etc. These studies along with others revealed prominent role of innate immune response to COVID-19 infection [[8], [9], [10], [11]]. The two important themes are that COVID-19 consistently results in the upregulation of chemokines, and increased levels of proinflammatory cytokines such as IL-1β, IL-2, and IL-6, thus leads to tissue damage. The second is that SARS-CoV-2 triggered a robust IFN response, marked by the upregulation of several ISGs in the lungs.

Like other viruses, novel SARS-CoV-2 also utilizes host machinery for their growth and survival during infection. Systematic investigation of virus–host protein–protein interactions (PPIs) offers an effective way toward elucidating the mechanisms of viral infection and drug repurposing [12]. It has also been reported that SARS-CoV-2 can establish a higher rate of infectivity if it acquires combinatorial mutations at the S-protein/ACE2 interfacial residues [13]. Several studies have shown that targeting cellular antiviral targets could be considered as a novel strategy for the development of effective treatments for viral infections, including SARS-CoV [14], MERS-CoV [14], Ebola virus [15], Zika virus [16], and more recently in SARS-CoV-2 [17].

There are no vaccines or drugs approved for the novel COVID-19 infection yet, but more than 80 clinical trials have been launched to test coronavirus treatments. The current COVID-19 pandemic leads the scientific community to focus seriously on drug repurposing to tackle the COVID-19 infection [[18], [19], [20], [21], [22]]. Toward this goal, in this study, we have generated a human-SARS-CoV-2 interactome based on recently published RNA-Seq analysis of human adenocarcinomic alveolar basal epithelial (A549) cells infected with SARS-CoV-2, and identified disease-related functional genes that will provide the insights into the patho-mechanisms of COVID-19. Moreover, drug-protein interactions of these hub proteins, investigated in this study will leads to the identification of potential candidate drugs for repurposing.

Materials and methods

Data collection

We have downloaded the differential gene expression data of human alveolar adenocarcinoma cells infected with COVID-19 from Melo et.al study [23]. We have also downloaded all the possible drug-gene interactions data from Drug-bank database and Comparative Toxicogenomic Database (CTD), along with drug-disease interactions as well [24,25].

Identification of significantly dysregulated genes

Melo et.al [26] previously described 120 differentially expressed genes (DEGs) based on p-value cut-off (p-val < 0.05). Identification of significantly dysregulated genes on the basis of P-value is not quite precise. The amount of change in the expression of genes (log2fold change) should also be considered. We reanalysed the data with q-value cut-off (FDR < 0.05), as well as applied log2fold change cut-off (log2fold ≥ 2) for identifying significantly dysregulated genes from the data of 120 genes.

Gene ontology and pathway enrichment analysis

For comprehensively analyzing the biological function of DEGs, we used database for annotation, visualization and integrated discovery (DAVID) version 6.7 [27]. This database uses the GO and the Kyoto encyclopedia of genes and genomes (KEGG) analysis for analyzing the DEGs. The gene ontology analysis included the annotation at biological level, cellular level and molecular level. The pathways and functions with FDR < 0.05 were considered significant.

Protein-protein interaction study

Using the DEGs list of Melo et al. [26], protein-protein interaction network was constructed using STRING plugin of cytoscape tool [28]. STRING plugin uses co-expression, text-mining, gene fusion, neighbourhood, and experimental data for preparing PPI network. Cluego tool was used for identifying the biological roles of DEGs and DAVID tool was used for gene ontology (GO) analysis.

Hub genes identification

In this PPI network, each node indicates a protein and an edge indicates an interaction between proteins. We have also calculated the network's topological parameters such as degree centrality (k), betweenness centrality (C ), closeness connectivity and eccentricity using network analyzer plugin of the cytoscape tool. Hub genes are the highly connected genes having high correlation with other genes in the network. Any changes in the expression of hub genes have the potential to influence the major part of the network. Genes having high connectivity can transfer information rapidly in the network as compare to other genes. Top genes which have the higher degree of connectivity (k) and betweenness centrality value were considered as hub genes.

Degree centrality (k) is value assign to each node purely based on number of links held by it. Degree indicates the number of interactions held by a node with other nodes in the network and helps in measuring the node significance in controlling the network.where, K is the node set containing all the neighbours of node u, and w(u,v) is the edge weight connecting node u with node v.

Betweenness centrality (C ) measures how many times a node falls on the shortest path with other neighbouring nodes. It characterizes a nodes ability to control the signal processing and information flow in the network.where p(k,u,f) is the number of interactions from k to f that passes through u, and p(k,f) represents total number of shortest interactions between node k and f.

Biological process enrichment analysis and GO analysis of the hub genes were done using Cluego and DAVID tool respectively.

Drug-gene interaction analysis

Drug–target interaction information for the hub genes was collected from the DrugBank database [24], and Comparative Toxicogenomics Database (CTD) [25]. The predicted drugs for hub genes through the gene-drug interaction databases were used for constructing drug-protein network using STITCH database [29]. STITCH is a database known to predict functional and physical interaction between chemicals/drugs and genes. The interactions in STITCH database is derived from five main sources, mainly by automated text-mining, high-throughput lab experiments, co-expression interaction data, interaction prediction by genomic context and by previous knowledges from databases. For each interaction between drug and gene, a combined score was calculated by STITCH database. The combined score was calculated by combining the probabilities of interaction from different evidence channels and corrected probability of randomly observing an interaction. The drug-gene interactions having network score more then 0.9 was consider significant and positively hit.

Results

Generation of human- SARS-CoV-2 Interactome by using transcriptome study

It is likely that the outcome of SARS-CoV-2 infection can largely be determined by the interaction between the host proteins and virus proteins. To build the human-SARS-CoV-2 interactome, we have used the transcriptome data of human alveolar adenocarcinoma (A549) cells infected with COVID-19 from Melo et al. [23]. Melo et al. identified 120 dysregulated genes with p-adj < 0.05, in which majority of genes were getting upregulated. Gene Ontology (GO) analysis of dysregulated genes using DAVID [30] revealed that the dysregulated genes were enriched in regulation of defence response to virus, regulation of protein export from nucleus, response to interferons, and regulation of response to cytokine stimulus (Fig. 1A). KEGG pathway enrichment analysis reveals the role of dysregulated genes in chemokine signaling pathways, complement and coagulation pathways and RIG-I-like receptor signaling pathways (Fig. 1B).

Fig. 1

Enrichment analysis of significantly upregulated genes. (A) GO pathway showed that upregulated DEGs were particularly enriched in which Biological process, Molecular functions, and Cell component. (B) The most significantly enriched KEGG pathway of the upregulated DEGs. (C) Gene network of 16 upregulated genes derived from GeneMANIA along with functional enrichment.

We reanalysed the transcriptomic gene expression signature based on log2fold change >2 and adj-p value < 0.05. We have identified 16 significantly dysregulated genes and of these genes, IFIT1, IFITM1, IRF7, ISG15, MX1, and OAS2 were highly upregulated. The GeneMANIA webserver was used to predict interactions between these DEGs genes in the network using the GO term “biological process” and source organism Homo sapiens as additional parameters. The gene enrichment of theses DEGs highlighted GO terms including response to virus, ribonucleotide binding and IFN-related signaling pathway (Fig. 1C). Overall, the analysis demonstrated that the upregulated genes are mainly linked to the host response to SARS-CoV-2 infection, type I interferon signaling and the cytokine-mediated signaling pathway.

We used the PPI network analysis to construct the interactome of the DEGs for COVID-19. For this, we used the STRING plugin of cytoscape tool [28], and obtained a PPI network with 117 nodes and 905 edges from the 120 dysregulated genes (Fig. 2A). The proteins are ranked based on their degree connectivity and betweenness centrality scores (Supporting information Table 1). Among these, STAT1, IRF7, IFIH1, MX1, ISG15, IFIT3, OAS2, and DDX58 has high degree and betweenness centrality values.

Fig. 2

Human-SARS-CoV-2 interactome. (A) The human PPI network analysis of the 120 upregulated genes, visualized by Cytoscape tool. Nodes with bigger size are considered as hubs in this network and are represented in diamond shape with red colour. (B) Functional enrichment analysis of the upregulated 120 genes that characterized the host response of COVID-19.

To obtain a more in-depth understanding of the interactome, GO function and KEGG pathway analysis were applied using DAVID (Fig. 2B and Supporting information Table 2). GO analysis results showed that in the biological process, the interacting genes were mainly enrichment in the regulation of defence response to virus, innate immune response, inflammatory response, and also played an important role in complement activation, and response to nutrient and extracellular stimulus (Fig. 2B). The cellular components are significantly located in the extracellular region and membrane fraction, etc. Molecular functions were mainly enriched in chemokines or cytokine activity, RNA binding, transferase activity, and calcium ion binding. KEGG pathway enrichment analysis revealed the role of interacting genes in complement and coagulation pathways, RIG-1-like receptor signaling pathways, and chemokine signaling pathways (Supporting information Table 2).

Overall, the network analysis indicates that SARS-CoV-2 targets the proteins in the IFN signaling pathway to evade the immune system. This highlights the key role of the IFN-mediated antiviral responses.

Network topology-based functional analyses reveal that SARS-CoV-2 targets core signaling pathways of the host network

Several virus-human interactome have shown that viral proteins significantly interacts with hub proteins of a network, and thus the topological features of a network can be used to predict viral targets [[31], [32], [33], [34]]. In this regard, two topological features, degree (number of connections) and betweenness (the fraction of all shortest paths that include a node within a network), were calculated to identify candidate hub nodes using the Cytoscape's Network analyzer tool.

Collectively, we have identified top 15 nodes with high degree of connectivity and Betweenness value and was subsequently considered as hub genes in the network (Table 1 ).

Table 1

Top 15 proteins based on the topological properties in PPI network of 120 dysregulated genes, and arranged on the basis of the decreasing degree of connectivity value.

Hub genes	Degree	Betweenness centrality
STAT1	55	0.247
IRF7	48	0.041
IFIH1	45	0.044
MX1	45	0.013
ISG15	44	0.051
IFIT3	44	0.028
OAS2	44	0.012
DDX58	43	0.015
IRF9	42	0.019
IFIT1	42	0.007
OAS1	42	0.006
OAS3	42	0.006
DDX60	41	0.027
OASL	41	0.006
IFIT2	41	0.003

Biological process enrichment analysis of 15 hub genes using GeneMania webserver and GO analysis revealed their role in defence mechanism to viral response, cellular response to type I interferon, regulation of viral genome replication, double-stranded RNA binding, type I interferon signaling pathway, and regulation of innate immune response (Fig. 3 and Supporting information Table 3).

Fig. 3

The PPIs network and gene enrichment analysis of 15 highly interactive Interferon stimulating genes (ISGs) along with functional enrichment.

These identified hub genes belong to the family of Interferon stimulating genes (ISGs). ISGs including IFIT and IFITM, ISG15, IFIH1, MX1, IRF7, OAS 1-3 and STAT1are known to potentiate IFN signaling and thus exert antiviral activity. These ISGs fall into two categories, one that is direct effectors of the innate immune response including: MX1, IFITM3, ISG15, IFIH1, and IRF-7, and the other includes the induction of viral RNA sensors such as DDX58 and the OAS1-3 genes. In addition to antiviral activity, ISGs may exert diverse functions including RNA and nucleotide binding, IFN regulation and inflammation regulation.

MX1, OAS1-3, IFIH1, IFITM1-3, ISG15 etc. are highly expressed in respiratory airways [35,36] and are associated with viral entry-associated ACE2 gene, as shown in single-cell RNA-sequencing data from various tissues from human [37,38]. Given the high expression of the viral entry-associated genes, it is reasonable that these nasal epithelial cells are trained to express these immune-associated genes to decrease viral susceptibility.

Further, we used Enrichr [39], a web-based server for gene-set enrichment analysis and provides different summaries of collective functions of gene lists. Interestingly, disease enrichment analysis demonstrated that the signature was highly associated with other viral diseases like West nile encephalitis, Tick-borne encephalitis, Dengue fever, Chikungunya, and SARS. Also, these hub genes are involved in the mammalian phenotype including increased susceptibility to viral infection induced morbidity/mortality, decreased interferon secretion, increased IgG level, abnormal T cell activation, and decreased interleukin-6 secretion.

Overall, these results clearly showed that the signatures are largely involved in innate immune response following SARS-CoV-2 infection (Supporting information Table 4).

Network based drug repurposing suggests toll like receptor 3 (TLR3) agonist as potential drugs to target COVID-19

Further using Drugbank database [24] and Comparative Toxicogenomics Database (CTD) [25], we identified the possible drugs which are known to have possible interaction with the hub genes (Supporting information Table 5). In total, we have identified 185 drug molecules which are known to interact with the hub genes (Supporting information Fig. 1). Apart from protein-drug interaction, we also looked for drug-disease interaction, i.e. which drug is used in which disease.

Based on the drug-gene interaction results, we used STITCH database [29] for final categorization of drug-gene interaction network based on interaction score which is ≥0.9. Interestingly, we found five proteins MX1, OAS1, STAT1, DDX58, and ISG15 were highly correlated to rare disease phenotype as well as to SARS in Enrichr database. These five proteins were also the most common influential proteins according to GO analysis. The STITCH drug-protein interaction network analysis of these five proteins has shown significant interaction with drugs/compound (Fig. 4 ). As shown in Fig. 4, MX1 gene interacts with Mitomycin-C, an antitumor and imiquimod, an immune response modifier. IFIH1 interacts with Polyinosinic:polycytidylic acid (poly I:C), which is an immunostimulant. OAS1 showed significant interaction with S-carbamidomethylcysteine (Cysteine-S-acetamide) and with MgATP. STAT1 interacts with Vanadium oxide and MgATP whereas, DDX58 and ISG15 interact with MgATP. Some of these identified drugs are used for induction of IFNs and thus play a key role in the body defence against SARS-CoV-2s infection.

Fig. 4

Drug-gene interaction network of five important ISGs and their related drugs using STITCH database. Green colour interactions between drug and genes were significant interaction with combined score > 0.9. The network shows poly (I:C), Imiquimod, Mitomycin C and other related drugs to ISGs network.

Discussion

Here, we generated a comprehensive human-SARS-CoV-2 interactome from transcriptome studies of a lung cell infected with CoVID-19. The PPI network analysis indicates that the pathways are enriched in host response to virus infection, type I interferons signaling, and cytokine activation. Network topology analyses identified 15 high-value targets of SARS CoV-2, which belongs to a subset of canonical ISGs. These ISGs are largely involved in regulation of defence response to virus, innate immune response, inflammatory response, and RNA binding [40,41].

An interferon-inducible protein, MX Dynamin Like GTPase (MX1) is associated with influenza and viral encephalitis infection [42,43]⁠. Also notable is the role of the Interferon Induced Protein with Tetratricopeptide Repeats 1 (IFIT1) and DExD/H-Box Helicase 58 (DDX58) as an antiviral activity. In Hepatitis E virus infection, polymerase binds to IFIT1 protecting the viral RNA from translation inhibition mediated by IFIT1 that boosts the interferon response in murine macrophage-like cells [44]. Recently, the significant upregulation of IFIT1-3, and DDX58 gene expression under COVID-19 viral infection has been reported [26]⁠. In response to viral infections, several genes of host such as OAS1-3, IRF7, IRF9, STAT1 and IFIH1 are highly expressed and are highly correlated with host response to viral infections [[45], [46], [47]].

Studies show that CoVs are equipped with strategies to antagonize the IFN signaling pathway that facilitates the virus to escape host immune response. SARS-CoV escapes the host IFN signaling as its ORF6 protein blocks the expression of STAT1-activated genes [48]. Similarly, in MERS-CoV ORF4b inhibits IRF3 and IRF7 to antagonize the antiviral IFN-β response [49]. Additionally, papain-like proteases (PLPs) are expressed in both SARS-CoV and MERS-CoV that enables to delay the host immune response. Coronaviruses engages in interactions with IFN stimulated gene 15 (ISG15) and antagonizing the IFN-mediated antiviral response [50,51]. The antiviral activity of ISG15 has been shown in several viruses including human cytomegalovirus [52], HIV [53], West Nile virus [54], swine fever virus [55], MERS-CoV [56], and SARS-CoV [57].

IFNs play a key role in the body defence against viral infections, and in this regard, we here showed some candidate drugs for repurposing. One of these drugs is poly (I:C) (polyinosinic:polycytidylic acid), a synthetic double-stranded RNA immune-stimulant, which is used as adjuvant in vaccine production [58]. It is agonist for toll like receptor 3 (TLR3) which induces the expression of IFNs. Many studies demonstrated that the TLR3 agonists, poly ICLC and poly (I:C) increases the production of IFN-α, -β, and -γ, which inhibited CoV replication and minimized the inhibitory effects of CoV on IFN signaling pathways [[59], [60], [61]]. Interestingly, chloroquine, a recently proposed drugs for COVID-19 [62], inhibits poly (I:C)-mediated IFN-β induction [63]. Therefore, TLR3 agonists can be a considered as potential drugs for repurposing in COVID-19. In addition to TLR-3 agonists TLR7 agonists such as imiquimod [64], can induce IFN production in the human body as well. Imiquimod is a strong inducer of IFN-α and several proinflammatory mediators, including TNF-α, IL-12, and chemokines [65].

The other drug predicted to bind MX1 and OAS, Mitomycin C is a cancer drug that is used in the treatment of bladder, colon, and breast cancers. It functions as an alkylating agent that causes cross-linking of DNA and inhibits RNA as well as protein synthesis [66,67]. At high concentration, Mitomycin C inhibits the replication of influenza virus by blocking the RNA synthesis [68]. Mitomycin C has also been shown to inhibit B cell, T cell, and macrophage proliferation in vitro and impair antigen presentation, as well as the secretion of IFN-γ, TNF-α, and IL-2. Vanadium oxide, an activator of STAT-1, involved in immune-regulating mechanisms, including immune suppression and inflammation downregulation by stimulating and activating B/T cells.3 [69,70]

Conclusion

In summary, our integrative interactome and network topology analyses showed that SARS-CoV-2 induced a strong IFN response, marked by the increased expression of several ISGs, including MX1, OAS1-3, IFIH1, ISG15, IRFs, and IFITMs etc. These ISGs exert antiviral activity and could protect the host cells from the infection. This protective action of ISGs might account for the lesser percentage of severe cases and the lower fatality rate in COVID-19. However, the extent of protection or damage to the host cell depends on stage of infection, types of cells, SARS clade [[71], [72], [73]] and other factors like co-infection, age, and co-morbidities. Recently, Zou et al. [74] reported high SARS-CoV-2 loads very early during infection, suggesting that the virus may have developed arsenals that is able to delay the IFN response by inhibiting innate immune signaling. Thus, IFN induction in the incubation period and at the very early stages of the infection could be the key to prevent COVID-19 associated mortalities. Administration of interferon-inducing agents, such as Poly (I:C) and Imiquimod could reduce the mortality of SARS at the very early stages of the disease (Fig. 5 ). On the other hand, at the later stages of the disease, the balance of the immune system becomes impaired, leading to inflammatory over-reactions, cytokine storm, and possible autoimmune responses. In such circumstances, administration of mitomycin C and other immunosuppressants to reduce possible lung inflammations may be needed (Fig. 5).

Fig. 5

Innate immune response as a protector to SARS-CoV-2. In the early stage of infection, SARS-CoV-2 inhibits host type I IFN antiviral immune defence. We propose to counteract this by utilizing toll like receptors (TLRs) agonists, poly (I:C) and Imiquimod to induce the production of IFNs and ISGs, which will reinstate the impaired immune responses. In the late stage of infection, there is a proinflammatory cytokine storm, which can be theoretically targeted by immunosuppressors, like Mitomycin C, and Vanadium oxide.

The following are the supplementary data related to this article.

Table S1

The human PPI network analysis of the 120 dysregulated genes.

Table S2

Functional enrichment analysis of 120 dysregulated genes using DAVID tool and with p-value ≤ 0.05 filter.

Table S3

Functional enrichment analysis of 15 hub genes.

Table S4

Gene set enrichment analysis for Molecular phenotype predictions in 15 Hub genes using Enrichr.

Table S5

Drug-gene interactions for 15 hub genes.

Supplementary Fig. 1

Network analysis of drug-genes of 15 hub genes along with the disease-drug interactions. Green triangles represent hub genes, Red square represents drugs and Blue circles represent diseases.

CRediT authorship contribution statement

Kartikay Prasad:Conceptualization, Methodology, Data curation, Software.Fatima Khatoon:Data curation.Summya Rashid:Writing - original draft.Nemat Ali:Visualization, Investigation.Abdullah F. AlAsmari:Supervision.Mohammad Z. Ahmed:Software, Validation.Ali S. Alqahtani:Software, Validation.Mohammed S. Alqahtani:Writing - review & editing.Vijay Kumar:Conceptualization, Supervision, Writing - review & editing.