Matrix metallopeptidase 9 as a host protein target of chloroquine and melatonin for immunoregulation in COVID-19: A network-based meta-analysis.

AIMS: The molecular pathogenesis of COVID-19 is similar to other coronavirus (CoV) infections viz. severe acute respiratory syndrome (SARS) in human. Due to scarcity of the suitable treatment strategy, the present study was undertaken to explore host protein(s) targeted by potent repurposed drug(s) in COVID-19.
MATERIALS AND METHODS: The differentially expressed genes (DEGs) were identified from microarray data repository of SARS-CoV patient blood. The repurposed drugs for COVID-19 were selected from available literature. Using DEGs and drugs, the protein-protein interaction (PPI) and chemo-protein interaction (CPI) networks were constructed and combined to develop an interactome model of PPI-CPI network. The top-ranked sub-network with its hub-bottleneck nodes were evaluated with their functional annotations. KEY
FINDINGS: A total of 120 DEGs and 65 drugs were identified. The PPI-CPI network (118 nodes and 293 edges) exhibited a top-ranked sub-network (35 nodes and 174 connectivities) with 12 hub-bottleneck nodes having two drugs chloroquine and melatonin in association with 10 proteins corresponding to six upregulated and four downregulated genes. Two drugs interacted directly with the hub-bottleneck node i.e. matrix metallopeptidase 9 (MMP9), a host protein corresponding to its upregulated gene. MMP9 showed functional annotations associated with neutrophil mediated immunoinflammation. Moreover, literature survey revealed that angiotensin converting enzyme 2, a membrane receptor of SARS-CoV-2 virus, might have functional cooperativity with MMP9 and a possible interaction with both drugs. SIGNIFICANCE: The present study reveals that between chloroquine and melatonin, melatonin appears to be more promising repurposed drug against MMP9 for better immunocompromisation in COVID-19.

Introduction

The pandemic ‘coronavirus disease 2019’ (COVID-19) is a severe respiratory illness caused by human coronavirus (HCoV), also known as SARS (severe acute respiratory syndrome)-CoV-2 or novel coronavirus which contains a single-stranded RNA genome [[1], [2], [3]]. The virus can primarily produce influenza-like symptoms viz. fever and myalgia along with pneumonia-like symptoms viz. dry cough and shortness of breath; however in extreme condition acute respiratory distress syndrome (ARDS), multiorgan failure and death may result [2,[4], [5], [6], [7]]. The transmission of SARS-CoV-2 from human to human occurs mostly by physical contacts, nasal droplets, uncooked foods and excreta of conciliator animals [1,2].

The SARS-CoV-2 phylogenetically belongs to other HCoVs under the genera of beta-coronaviruses viz. epidemic SARS-CoV and MERS (Middle East respiratory syndrome)-CoV. The whole genome of SARS-CoV-2 has 79% nucleotide sequence homology to that of SARS-CoV. The surface envelope spike glycoprotein-S is the major antigen of both SARS-CoV and SARS-CoV-2 having 75% amino acid similarity. The glycoprotein-S binds to its receptor, a plasma membrane bound protein called angiotensin converting enzyme subtype 2 (ACE2) of human, allowing endocytosis of the complexes and entry of viruses into the host cells [1,2,8]. The host defense mechanism during HCoV infection gets compromised, as the immune defense responses against the interferons occur late which ensures enough time for replication and survival of the viruses [9]. Clinically, interferon supplementation in COVID-19 patients has been reported to reduce the virus load and inflammatory insult [10].

Extensive review analyses using non-randomized clinical studies [6,[11], [12], [13], [14], [15]] and in vitro/culture studies [[16], [17], [18]] reveal that several approved drugs including anti-viral (remdesevir) [13,14,16,17], anti-retroviral (ritonavir, darunavir, lopinavir) [13,15], anti-malarial (chloroquine, hydroxychloroquine) [13,14,16,17,19], anti-protozoal (nitazoxanide, ivermectin) [14,18] and immunosuppressive (tocilizumab) [13] agents may be the choice of treatment for COVID-19. Monoclonal antibody against the spike protein of SARS-CoV-2 has been raised in view of developing protection against COVID-19 [20]. Moreover, the network analysis based on the human coronavirus genome and host protein interactome has recently proposed that some repurposed drugs including anti-inflammatory agent (melatonin), non-steroidal selective estrogen receptor modulator (toremifene), angiotensin receptor blockers (e.g. irbesartan), immunosuppressants (e.g. sirolimus), anti-neoplastic drugs (e.g. mercaptopurine) and a natural plant product (emodin) having antiviral activities may be the likely candidates for treatment of COVID-19 [21]. However, all these pharmacologic agents suffer from certain limitations in their uses.

Due to the paucity of information about the molecular pathogenesis of COVID-19 and its best suitable treatment, the present study has been executed with potential systems biology approach using integrative protein-protein and chemo-protein network analyses to find out (a) the most effective host protein target(s) in human and (b) the putative drug(s) against the target(s), to establish a functional link between drug(s) vs. targeted host protein(s) for better understanding of treatment strategy in COVID-19.

Materials and methods

The systematic and stringent methodology with inclusion and exclusion criteria applied in the present study is given in the flow diagram (Fig. 1 ).

Fig. 1

Flowchart of the systematic and stringent methodology applied and the results found in the mechanistic systems biology analysis to identify potential drugs against their targeted biomolecule(s) in COVID-19.

Raw data acquisition and processing to identify differentially expressed genes (DEGs)

The microarray dataset (GSE1739) of gene expression profiles in the blood of 10 SARS-CoV patients along with that of four healthy individuals were collected from the NCBI Gene Expression Omnibus (GEO) database. Diagnostic profiles and experimental assay with the including and excluding selection criteria for human subjects had been described in the original publication [22]. The basic packages of the Bioconductor project [23] were used for retrieval, background correction (logarithmic transformation) and quantile normalization of the data in R language and environment [24]. A 95% confidence interval was applied to select the data for further analysis. Further, the R package limma [25] was implemented to analyze the expressions of genes in patients and healthy individuals. The lmFit function was used to fit the linear model using the least square method followed by another function eBayes (a hierarchical Bayesian model) for computation of moderated t-statistics on the linear model to find differential gene expression. The moderated t-statistics uses an empirical Bayesian shrinkage estimator to reduce the variance dependency for a specific gene on its mean expression values [26]. The differential expression of genes was computed using ‘false discovery rate’ (FDR) based on Benjamini and Hochberg method [27] with the criteria of adjusted p-value < 0.05 and |log2(FC)| > 1 (FC representing ‘fold change’) for each comparison.

In silico modelling of physical ‘protein-protein interaction’ (PPI) and ‘chemo-protein interaction’ (CPI) networks

The DEGs found in the present study, had been incorporated on the STRING 11.0 database [28] with the settings for active interaction from all enable sources (‘Textmining’, ‘Experiments’, ‘Databases’, ‘Co-expression’, ‘Neighborhood’, ‘Gene Fusion’, ‘Co-occurrence’). The reliable interaction strength was detected by applying parameter customized confidence score of 0.600 as threshold to filter out spurious interactions. The study included only the moderate and high probable interactions between the protein partners of differentially expressed genes to construct an integrative physical ‘protein-protein interaction’ (PPI) network model for the characterization of molecular interactions involved in SARS-CoV infection.

Then the systematic review had been performed to acquire recent relevant literatures associated with COVID-19 therapy in PubMed database () using the following keywords: ‘SARS’,‘SARS-CoV’, ‘COVID-19’, ‘treatment’, ‘therapy’, ‘therapeutics’ and/or ‘Drug’ respectively. The search results were filtered with ‘Publication date 1 year’ for the latest updates on COVID-19 therapies to select drug compounds from the latest literatures [18,21,29,30]. The drugs (targeted to COVID-19/SARS-CoV-2 treatment) selected and DEGs (identified from the peripheral blood samples of SARS-CoV patients) found in the present study were submitted to STITCH 5.0 database [31] with the setting of active interactions from all enable sources (‘Textmining’, ‘Experiments’, ‘Databases’, ‘Co-expression’, ‘Neighborhood’, ‘Gene Fusion’, ‘Co-occurrence’). The known and predicted binding interactions between genes/proteins and small chemical compounds/drugs were detected by applying parameters including customized confidence score of 0.600, ‘interactor/query protein only’ and ‘network depth equal to 2’ to construct a ‘chemo-protein interaction’ (CPI) network that included moderate and strong target (protein)-drugs interactions for COVID-19.

Furthermore the CPI network model (developed in STITCH) and PPI network model (developed in STRING) were imported in Cytoscape 3.7.2 software [32] and merged by applying the ‘union function’ of the Cytoscape core plugin ‘Merge Networks’ based on attribute values of nodes and edges to get interactome model of merged physical ‘protein-protein and chemo-protein’ (PPI-CPI) network to identify both the potential drug targets (proteins) and potential drug candidates for COVID-19.

Topology analysis of interactome model of PPI-CPI network to find sub-networks and hub-bottleneck nodes

Molecular Complex Detection (MCODE) [33] (a Cytoscape plugin available at ) module using cut-off parameters (scoring and cluster finding) and the MCODE score ≥ 4 was utilized to screen the significant/major sub-network of PPI-CPI interactome for COVID-19. The scoring cut-off parameters were ‘loop included’ (inclusion of loops i.e. self-edges in neighborhood density calculation) and ‘degree cut-off 2’ (number of connections necessary for a node to be scored). The cluster finding cut-off parameters were ‘node score cut-off 0.2’ (cluster size for selection of new members), hairfall included (deletion of all single connected nodes from clusters), fluff included (allowance of cluster expansion with one neighbor shells), node density cut-off 0.1 (controls the neighbor inclusion criteria during ‘fluffing’), K-core cut-off 2 (filters out clusters by eliminating small inter-connected sub-cluster) and maximum depth of network cut-off 100 (limiting the distance from the seed node within which the cluster member can be searched).

The chemo-protein components of overall PPI-CPI interactome were analyzed in CentiScaPe 2.2 [34] plugin of Cytoscape using key network centrality parameters like node degree (the number of neighboring nodes to which the node of interest is directly connected), shortest path betweenness (the number of information streams passing through a given node) and stress (the extent to which a node can hold network communications) to identify the most potential nodes in the PPI-CPI network. Most relevant nodes were selected using the numerical threshold values according to the mean centrality values for each of the topological centrality properties: ‘mean node degree’: 5.02; ‘mean betweenness’: 277.21; ‘mean stress’: 934.0. Nodes (protein/chemical) with high degree/connectivity centrality, high betweenness centrality, and high stress centrality (above mean threshold values) were considered as the hub-bottleneck nodes, which were finally taken into account to serve as the functionally essential nodes (genes/proteins/drugs) of this network model.

Functional enrichment analysis of the interactome model of the top-ranked sub-network

The list of genes of interest with their identifiers (ID) found in the top sub-network of PPI-CPI interactome was submitted to Enrichr web-tool platform () [35,36] for the enrichment (i.e. over representation of common annotated biological features) analysis of functional annotations implemented across a number of resources including ‘gene ontology’ (‘biological process’, ‘cellular component’, ‘molecular function’), ‘KEGG biological pathways’, ‘Jensen Disease’. The Enrichr analysis was performed using statistical parameters viz. p-value (Fisher exact test), q-value (adjusted p-value for false discovery rate), old p-value, adjusted old p-value, odd ratio, z-score and combined score (log(p-value) × z-value). The Enriched results (functional annotations) were ranked based on the levels of significance with p-values < 0.05 and corresponding combined scores followed by selection of terms with certain cut-off (top 10 terms).

Results

The results found in the present mechanistic systems biology analysis are systematically documented in the flow diagram with the findings of potential drugs against the targeted biomolecule(s) in COVID-19 (Fig. 1).

Identification of DEGs using microarray data of SARS-CoV patient blood

From in silico analysis of SARS microarray dataset we identified 120 differentially expressed genes among which 45 genes were upregulated and 75 genes were downregulated. Details of the differentially expressed genes including gene identifiers (ID), “false discovery rate” (FDR) adjusted p-values (<0.05) and log2(fold change) values (>1) are summarized in Fig. 2 .

Fig. 2

Volcano plot analysis to identify the differentially expressed genes (DEGs). The expressions of genes are evaluated by analysis of microarray data of the peripheral blood samples of SARS-CoV patients (n = 10) versus healthy controls (n = 4) collected from the data source (GSE1739) using Bayesian algorithm in limma Bioconductor package of bioinformatics tools in R language and environment. A: The volcano plot of the expressions of genes using the logarithmic values of fold changes (log2 fold change) in x-axis and ‘false discovery rate’ (FDR) adjusted p-values (−log10 adjusted p-value) in y-axis. Dotted lines parallel to x-axis and y-axis indicate the threshold values using FDR adjusted p-value < 0.05 and | log2 fold change | > 1 respectively to identify the upregulated and downregulated DEGs. The red, blue and black color dots indicate the upregulated DEGs, downregulated DEGs and non DEGs respectively. B: The gene IDs of upregulated (upper panel) and downregulated (lower panel) DEGs found in volcano plot. The position of dots of DEGs (45 upregulated and 75 downregulated) are same as those appeared in volcano plot. The scales of upregulated and downregulated DEGs are adjusted manually for proper presentation of the gene IDs. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Construction of interactome model of PPI-CPI network for COVID-19

The interactome model of PPI network was constructed using proteins corresponding to the respective DEGs of SARS-CoV patients in STRING webtool. A total of 72 protein nodes and 212 connections had been found in PPI network (data not shown). Literature survey provided a total of 65 potential drugs proposed for treatment of COVID-19 (Table 1 ) and these drugs were included in further analysis for selection of most potent drug(s) against putative protein target(s). The interactome model of CPI network was constructed using COVID-19 drug candidates and DEGs of SARS-CoV patients in STITCH webtool. Total 88 nodes (proteins and drugs) and 120 connections had been found in the CPI network (data not shown).

Table 1

Summary of the potential drug candidates selected from recent literature on COVID-19. Drugs are categorized on the basis of their mode of actions.

Category	Drug candidate(s)
Analgesic	Diperodon, Phenazopyridine, Tetrandrine
Anti-bacterial	Dihydrocelastryl diacetate, Monensin sodium, Oligomycin, Salinomycin sodium, Valinomycin
Anti-depressant	Desipramine
Anti-fungal	Antimycin A, Exalamide, Phenylmercuric acetate
Anti-helminthic	Pyrviniumpamoate, Ivermectin
Anti-histamine	Chloropyramine
Anti-hypertensive	Alprenolol, Berbamine, Carvedilol, Doxazosinmesylate, Irbesartan, Propranolol
Anti-infective	Cetylpyridinium chloride, Camphor
Anti-inflammatory	Colchicine, Emodin, Mesalazine
Anti-malarial	Chloroquine, Conessine, Hydroxychloroquine, Quinacrine
Anti-neoplastic	Dactinomycin, Hydroxychalcone, Lycorine, Mercaptopurine, Mycophenolate mofetil, Mycophenolic acid, Pristimerin, Toremifene
Anti-Parkinsonian	Harmine
Anti-protozoal	Nitazoxanide
Anti-psychotic	Promazine
Anti-viral	Acyclovir, Favipiravir, Ganciclovir, Lopinavir, Oseltamivir, Penciclovir, Remdesivir, Ribavirin, Ritonavir, Tilorone
Ca2+ channel blocker	Loperamide
Diuretic	Eplerenone
Estrogen steroid	Equilin
Hormone	Melatonin
IL-6 inhibitor	Tocilizumab
Immunosuppresants	Sirolimus
Muscle relaxant	Papaverine, Zoxazolamine
PAF inhibitor	Ticlopidine
Protein synthesis blocker	Cycloheximide, Emetine
Selective serotonin reuptake inhibitor	Paroxetine
Steroid hormone	Oxymetholone
Serine protease inhibitor	Nafamostat

The characterization of protein-drug interactions in COVID-19 had been executed through development of interactome model of PPI-CPI network in Cytoscape by merging the already identified PPI and CPI network and that was found to be composed of 118 nodes (proteins and drugs) interconnected by 293 interactions (Fig. 3A). The identification of nodes of gene products/proteins and drugs has been designated by the corresponding gene IDs and name of the drugs respectively in the present study. Accordingly, the respective gene IDs have been used as the node (protein) identifiers in the further descriptions.

Fig. 3

The interactome models of PPI-CPI network (A) and its top ranked sub-network (B) obtained in the Cytoscape software. The networks are developed using DEGs of the SARS patients and the promising COVID-19 repurposed drug candidates. A: The PPI-CPI network consists of 118 nodes of gene products/proteins and drugs corresponding to their respective gene IDs (35 upregulated and 43 downregulated DEGs) and name of drugs (40 in number). The 293 edges correspond to the functional connectivities between nodes. The upper panels of the images represent networks having continuous connections. The lower panels of the same designate discrete networks having only single connections. B: The top ranked sub-network was identified in the MCODE module of the Cytoscape considering MCODE score ≥ 4. The sub-network consists of 35 nodes of gene products/proteins and drugs of corresponding gene IDs (26 upregulated and 5 downregulated DEGs) and name of drugs (4 in number) respectively. The 174 edges correspond to the functional connectivities between nodes. A and B: The shape of the nodes in interactome models (A and B) depicted as circles represent the nodes associated with the top ranked sub-network. The rest of the nodes in PPI-CPI network (A) are manually altered to diamond shapes to identify the top ranked sub-network in the main network for better interpretation. The intrinsic attributes of nodes and edges are kept intact in interactome models (A and B). The border of the nodes is illustrated with 2 pts. grey color. The sizes of the nodes indicate their connectivities (higher the value, higher will be the size) adjusted by the ‘continuous mapping of node size’ in ranges between 25 and 60 pts. for the lowest and highest node size respectively. The color of the nodes is represented as cyan for the drugs. The color gradients of red:white:blue are denoted as the expression pattern (log2 fold change) of genes specific for gene products/proteins (nodes), whereby the gradients of the colors have been adjusted to the expression values over the ranges of 6.68:0:−6.68 from the ‘continuous mapping of node color’ in the node network style of the Cytoscape. The widths of the edges are based on the combined scores (0.600 to 1) obtained as interaction data in the STRING for PPI network and STITCH for CPI network which are adjusted by 0.5 to 5 pts. of ‘continuous mapping of edge width’ in the edge network style of the Cytoscape. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Identification of top-ranked sub-network and hub-bottleneck nodes within the interactome model of PPI-CPI network for COVID-19

From the interactome model of PPI-CPI network, the MCODE module in Cytoscape identified only one top-ranked sub-network (cluster having MCODE score 9.33) that comprised 35 nodes (proteins and drugs) with 174 interconnections (Fig. 3B) among which 31 nodes appeared to be protein molecules and four viz. chloroquine, melatonin, propranolol and quinacrine were found to be drug candidates. The sub-network (Fig. 3B) contributed the 59.38% major connections and 29.66% nodes in the overall interactome model of PPI-CPI network (Fig. 3A).

The top-ranked sub-network showed hub-bottleneck nodes consisting of 12 proteins and drugs with high topological centrality index values (high node connectivity, high betweenness, and high stresses, each with greater mean cut-off threshold values) that include 10 proteins corresponding to genes viz. CAMP, CCT2, ELANE, FOXO3, ITGAM, MMP9, MPO, SIRT1, SMAD4, STAT1 and two drug candidates, namely chloroquine and melatonin (Fig. 4A–C). The hub-bottleneck nodes were identified using the cut-off thresholds of node degree (value 5.02), node betweenness (value 277.21) and node stress (value 934.0) applied to the top-ranked sub-network (Fig. 4A–B).

Fig. 4

Pictorial summary of the topological properties and the centrality analyses for the top ranked sub-network to identify the hub-bottleneck nodes (A, B, C) using the CentiScaPe module of Cytoscape software. The identification of nodes of gene products/proteins and drugs are designated by corresponding gene IDs and name of the drugs. Graphical plots represent the dot plots of values of (A) node degree (x-axis) vs. node betweenness (y-axis), (B) node degree (x-axis) vs. node stress (y-axis) and (C) Venn diagram of high node degree/connectivity, high node betweenness and high node stress. Here, the term ‘high’ indicates higher than the mean cut-off thresholds for node degree/connectivity, betweenness and stress, which have been obtained from the CentiScaPe module of the Cytoscape software. Mean centrality values are presented as dotted lines in the graphs (A, B). The black round dots are hub-bottleneck protein nodes, the blue diamond shapes are non-hub-bottleneck nodes (proteins and drugs) and red boxes with black borders are the hub-bottleneck drug nodes in the (A, B) graphs. MMP9 in the red outlined boxes in the graph (A, B) represent the only target of chloroquine and melatonin among other hub-bottleneck nodes of the top-ranked sub-network. (C) Venn diagram indicates the common nodes that have the topological centrality indices viz. node degree/connectivity, betweenness and stress with the values higher than the mean cut-off respective threshold values obtained from the CentiScaPe module of the Cytoscape software. The upward and downward arrows indicate the expressions of upregulated and downregulated genes (gene IDs right to the arrows) corresponding to respective gene products/proteins in the hub-bottleneck nodes of COVID-19. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Among the 10 hub proteins (Fig. 4C) six proteins corresponding to respective upregulated genes (CAMP, ELANE, FOXO3, ITGAM, MMP9, MPO) and four proteins corresponding to respective downregulated genes (CCT2, SIRT1, SMAD4, STAT1) (as evident in their expression profile in Fig. 2), only matrix metallopeptidase 9 (MMP9) protein had been found to interact with both hub drug nodes, namely chloroquine and melatonin, in the top ranked sub-network (Fig. 3B). MMP9 had been found as the second largest hub node consisting of connectivities with 23 nodes (19 protein and four drug partners). Notably, the first largest hub node ELANE had connectivities with 26 nodes without having any interacting drug partner (Fig. 3B). Furthermore, the other hub node FOXO3 appeared to have connectivity with melatonin as the only interacting drug partner (Fig. 3B). However, there is no evidence in support of the involvement of gene product of FOXO3 with the respiratory viral infection. Therefore our present results clearly indicated that MMP9 could be the most potent host protein target of chloroquine and melatonin as interacting drugs in COVID-19.

Functional enrichment analysis of the interactome model of the top-ranked sub-network for COVID-19

Gene enrichment analysis study (Fig. 5 ) indicated that the gene products of the top-ranked sub-network were associated with biological processes including neutrophil related processes (neutrophil activation, neutrophil degranulation, neutrophil mediated immunity and immune responses), innate immune responses in mucosa, humoral immune responses mediated by antimicrobial peptides, antibacterial humoral responses, defense responses to bacteria with high statistical significances (Enrichr p-value < 0.05 and combined score). The association of the gene products of this sub-network were found to be statistically significant (Enrichr p-value < 0.05 and combined score) for the cellular components involving subcellular and granular secretary functions. Statistically significant (Enrichr p-value < 0.05 and combined score) molecular functions like nuclear hormone receptor bindings, iron and transition metal ion bindings, serine-type endopeptidase activities, endonuclease and ribonuclease activities and RNA polymerase-II core promoter proximal region sequence-specific DNA bindings were also associated with the gene products of this sub-network. In addition, KEGG pathway enrichment study included certain signaling pathways (FoxO, NOD-like receptor and prolactin), several types of microbial infections (exerted by Staphylococcus aureus, Mycobacterium tuberculosis, Leishmania donovani and Hepatitis-B virus) and cancer (pancreatic and acute myeloid leukemia) in association with the gene products of the sub-network. Jensen disease enrichment data indicated that the respiratory system disorders (lung diseases, bronchitis and common colds), immunoinflammation related complications (Wegener's granulomatosis, Mastitis, Arthritis), eosinophilia, brain edema, periodontal disease were highly enriched in association with the gene products of the top-ranked sub-network.

Fig. 5

The bubble plot to demonstrate the summary of the functional enrichment analysis of DEGs under the top ranked sub-network cluster of the PPI-CPI network in COVID-19 identified by the Enrichr program. The results include the top ten enriched terms corresponding to each of the Enrichr functional annotations viz. GO terms [Biological Process (GO-BP), Cellular Component (GO-CC), Molecular Function(GO-MF)], KEGG pathways and Jensen Diseases obtained from Enrichr analyses using 31 genes corresponding to the gene products/proteins nodes of top ranked sub-network of the PPI-CPI network in COVID-19. The abscissa represents the ‘Enriched gene ratio’. The sizes of bubbles indicate the ‘Enriched gene count’ adjusted with 0.5 to 4 pts. based upon the gene counts of respective functional annotations. The colors of the bubbles indicate the ‘Enrichr combined scores’ (log(p-value) × z-value) which are adjusted in a color (VIBGYOR) gradient, ranging from 0 to 7000 on the basis of values of the scores. The p-value < 0.05 of Fisher's exact test is considered for significant result. The black colored asterisks indicate the enriched functional annotations of the gene products/proteins including MMP9 that correspond to their DEGs denoted by numbers (1 to 13). 1–3: ARG1, AZU1, BPI, CAMP, CCT2, CEACAM6, CEACAM8, CHI3L1, CRISP3, DEFA1, DEFA4, ELANE, GNS, ITGAM, LCN2, LTF, MMP9, MPO, MS4A3, PGLYRP1, RNASE2, RNASE3, SLPI, TCN1; 4: CAMP, CRISP3, LTF, MMP9, PGLYRP1, TCN1; 5: CAMP, CEACAM8, CRISP3, ITGAM, LTF, MMP9, PGLYRP1, TCN1; 6–7: AZU1, ELANE, MMP9; 8: RNASE2, RNASE3; 9: MMP9, MPO, RNASE3; 10: MMP9, MPO; 11: ITGAM, STAT1, HP, CHI3L1, MPO, ELANE, MMP9; 12: CEBPE, ELANE, ITGAM, MMP9, MPO; 13: MMP9, SMAD4, STAT1.

Discussion

The present study developed a chemo-protein interactome network (Fig. 3A) on basis of differentially expressed genes in SARS-CoV infection (Fig. 2) and repurposed drugs for COVID-19 (Table 1). A sub-network (Fig. 3B) was identified that appeared to be involved in neutrophil activation and degranulation pathways (Fig. 5). The functionally important hub node MMP9, an immunogenic protease linked to the degradation of cellular matrix, was revealed as the target (Fig. 5) of two repurposed promising drug candidates (Figs. 3B and 4C) i.e. chloroquine (antimalarial drug) and melatonin (pineal hormone) in COVID-19. Here, the systematic analytical approach using data integration reveals the global systems-level relationship (drug vs. target vs. functional annotation) in COVID-19.

ARDS is the main cause of morbidities in COVID-19 and infections of other coronaviruses [[5], [6], [7],37,38]. Cytokine storm is a key mechanism of ARDS materialization leading to multiorgan failure and death [5,6,39,40]. Till date, the molecular pathogenesis of COVID-19 is unclear. Therefore, the similar mechanism of SARS-CoV/MERS-CoV infection can confer a lot to the molecular level of understanding in pathogenesis of COVID-19 [6]. The present study identified a sub-network (Fig. 3B) that includes the statistically significant pathways of both neutrophil activation and degranulation (Fig. 5). The sub-network comprised 29.66% nodes (proteins) and provided 59.38% connectivities of the overall network, which indicated the colossal importance of this sub-network in the whole network model (Fig. 3A–B). The global and local topological analyses provided ten potential protein molecules consisting of six proteins corresponding to respective upregulated genes (CAMP, ELANE, FOXO3, ITGAM, MMP9 and MPO) and four proteins corresponding to respective downregulated genes (CCT2, SIRT1, SMAD4 and STAT1). These ten nodes appeared to be the hub nodes of the PPI-CPI interactome network and were chosen for further study (Fig. 4C). All these molecules had been found to be highly enriched in blood after antimicrobial-induced neutrophil-mediated humoral and innate immunity responses in COVID-19 (Fig. 5).

In the present study, MMP9 was found as the functionally important hub node (protein) and target of the drugs in the central hub node of the PPI-CPI interactome network (Figs. 3B and 4C). Matrix metallopeptidases (MMPs), the zinc-containing and calcium-dependent proteolytic endopeptidases, are released from the intracellular stores and become active extracellularly [41]. MMPs cause deterioration of a number of extracellular matrix proteins that help in the extracellular matrix remodelling in various physiological and pathological processes including inflammation [[41], [42], [43], [44]]. MMPs are widely distributed in tissues and their expressions are controlled by cytokines, growth factors and hormones. MMP9 (83kD), a glycoprotein with type-IV collagenase activity, is found in neutrophils, lymphocytes and dendritic cells that are involved in angiogenesis and inflammatory cytokine generation [41,44]. Notably, the Enrichr functional analysis in the current study revealed MMP9 as to be associated with transition ion binding (molecular function under gene ontology) and diverse pathophysiological conditions viz. neutrophil activation and related immunity (biological process under gene ontology), transcriptional misregulation in cancer and hepatitis B infection (KEGG pathways), arthritis, brain edema and common cold infection (Jensen disease) in human (Fig. 5). The computational study using functional pathways enrichment analysis has also reported the bidirectional regulation of MMP9 gene expression in diverse viral infections and predicted MMP9 as the repurposing drug target [45].

MMP9 is constitutively expressed at low level under physiological condition. The MMP9 action depends on the balance between its pro-peptide activation by a serine protease (plasmin) and inhibition through formation of complexes with ‘tissue inhibitor of matrix metalloprotease 1’ (TIMP1). Therefore, the dysregulations of MMP9 systems (gene expressions, post-translational modifications, levels of activators vs. inhibitors) lead to the development of pathological conditions [46,47]. The low level of TIMP1 and persistent MMP9 activity in neutrophils which bring about airway remodelling and congestion in association with the acute inflammatory responses are implicated in development of both ARDS and acute asthma [[48], [49], [50], [51]]. Recently the gene expressions of TIMP1 in the bronchoalveolar lavage fluid of COVID-19 patients have been found to undergo upregulation in one victim and downregulation in another one as well as upregulation in the blood mononuclear cells of the other three victims [52]. In the present communication, the Enrichr functional analysis (cellular components database under gene ontology) exhibited the association of MMP9 with functional aspects of lumens releasing tertiary granules (Fig. 5). MMP9 is reported to relate with tertiary lumens for release of cytokines during neutrophil activation [53]. In monocytes and macrophages, MMP9 is over-activated by plasmin and reacts with ‘toll like receptor 9’ (TLR9) signaling to induce formation of tumour necrosis factor (TNF) that leads to the development of pro-inflammatory cytokine storm [54,55]. The involvement of MMP9 with serine protease activities was also noticed in Enrichr functional analysis with molecular function under gene ontology (Fig. 5). MMP9 deficiency is found to be protective against severe H1N1 influenza virus A infection in mice model [56]. Therefore, our findings and other reports support the view that the approach to MMP9 inhibition and/or alterations of its activators (plasmin) and inhibitors (TIMP1) may prevent the deadly cytokine storm and the life-risk of COVID-19 patients.

Here, both chloroquine and melatonin were identified as the repurposed drug candidates that appeared to be the interacting partners of MMP9 (Figs. 3B and 4). Chloroquine treatment decreases serum MMP9 level in systemic lupus erythematosus [57] and suppresses MMP9 activity and its mRNA expression in breast cancer [58]. Chloroquine can also lower the TNF-mediated neutrophil apoptosis, neutrophil degranulation and cytokine burst [[59], [60], [61], [62], [63]]. Chloroquine/hydroxychloroquine efficiently inhibits the human retrovirus activities in vitro [64]. Recently, a number of studies [17,65,66] and a clinical trial [67] support the use of chloroquine/hydroxychloroquine for the treatment of COVID-19. Conversely, certain clinical studies strongly disagree with its uses [[68], [69], [70]], as it develops cardiac failure in COVID-19 patients [68], indicating the clinical data about chloroquine/hydroxychloroquine are not conclusive [71]. Notably, hydroxychloroquine did not appear in the sub-network in the present study (Fig. 3B) and remained unattained for further consideration.

Melatonin is responsible for normal sleep and maintenance of the ‘biological clock’ in human [72]. It has anti-inflammatory properties [73,74], exerts its potential role in anti-viral mechanism [[74], [75], [76], [77], [78]] and has also been used in Ebola virus infection [78]. The pathophysiological observation in COVID-19 patients supports the chronobiological uses of melatonin in its treatment [79,80] and a therapeutic algorithm in this regard has also been proposed recently [81]. In addition, melatonin is reported as the promising adjuvant for COVID-19 treatment [21,82] and its deficiencies may enhance susceptibility of diabetic and hypertensive elderly patients in SARS-CoV-2 infection [80]. The first clinical trial of melatonin in a small cohort indicates its involvement in abating ARDS in COVID-19 patients [83]. Interestingly, melatonin can bind to the active site of MMP9 and inhibit the latter to arrest immunoinflammation [84]. Thus, the present network based meta-analysis justified the involvement of melatonin interaction with MMP9 in immunocompromised COVID-19.

The cell culture study (kidney cell line Vero E6) and network analysis reveal that chloroquine [85] and melatonin [21] respectively may have potential anti-viral role in targeting human cell membrane-bound ACE2 receptor, a zinc dependent carboxypeptidase, and other associated protein partners in COVID-19. Interestingly, separate studies report that MMP9 level increases in ACE2-knockout mice model [86] and ACE inhibitors (lisinopril and imidapril) also target MMP9 along with ACE2 [87]. Therefore, both metallopeptidases i.e. ACE2 and MMP9 may have cooperativities for immunoinflammation in COVID-19 and can be potential targets of both chloroquine and melatonin therapies.

The glycoprotein-S of surface spikes of corona viruses including SARS-CoV-2 recognize the membrane-bound ACE2 receptors in airways and lungs for entry into the body. ACE2 along with transmembrane serine protease (TMPRSS2) triggers conformational changes in the glycoprotein-S and releases S-fragments that enable the fusion of viral envelopes with the host cell membrane followed by internalization of viruses through the formation of endosomes. The viral RNA is processed within the host cells and replica of viruses are released out of the cells. The viral infected cells present antigenic peptides through major histocompatibility complexes and elicit humoral and cellular immunities along with cytokine storm [6,[88], [89], [90], [91], [92], [93]]. The ectodomain of membrane-bound ACE2 is shed into the extracellular fluid as its soluble form [88,[94], [95], [96]]. In vitro studies also indicate that the soluble ACE2 can bind to SARS-CoV/SARS-CoV-2 and is supposed to limit the availability of viruses to interact with the membrane-bound ACE2 receptors thus minimizing the chances of viral load [88,92,94]. Notably, ACE2 converts angiotensin-II to angiotensin(1-7), a peptide which causes modification of inflammatory processes [97]. Angiotensin(1-7) also has a vasodilatory effect [98] and enhances the heart rate acting on the caudal ventrolateral medulla in the brain [99], which may contribute to the development of hypotension and tachycardia as found in SARS/COVID-19 patients [100]. Moreover, the interaction of SARS/SARS-CoV-2 with the ACE2 receptors of the support cells in olfactory mucosa also results in inflammation-induced anosmia [101]. Further experimental studies are needed to elaborate the regulatory role of ACE2 receptor and MMP9 in therapeutic management of COVID-19.

Conclusion

The present study using network-based systems biology approach clearly demonstrates that (a) the upregulation of MMP9 gene and status of MMP9 protein with its activity/activity regulators, as potentiators of cytokine storms, may be the key point of etiopathophysiology of COVID-19 and (b) chloroquine and melatonin may be the drugs that can target MMP9 to reduce the immunoinflammatory cascades associated with ARDS symptoms in COVID-19 patients. Consistent with the present findings, it would also be worthy to mention that the therapeutic actions of both chloroquine and melatonin may have converged onto the activities of cell membrane ACE2 receptors of SARS-CoV-2. Recently the use of chloroquine does not appear attractive as it produces cytotoxicity at doses used in the clinical trials. In this respect melatonin could be a better promising safe drug that needs more clinical trials in large cohorts to ensure its efficacy in treatment of COVID-19 patients.

Declaration of competing interest

The authors declare that there is no conflict of interest.