Introduction

Coronavirus disease 2019 (COVID-19) has emerged as a public health concern in terms of mortality and economic loss. Originating from China, COVID-19 drastically spread throughout the world quickly and affected more than 50 different countries (Zhang et al. 2020). World Health Organization declared COVID-19 as “Pandemic” and has suggested preventive approaches like social isolation and intensive care of patients (WHO 2020). Recent in-silico and network pharmacology studies have reported repurposing of existing drugs and the alleged role of natural lead in the treatment of COVID-19 (Sinha et al. 2020a,b; Chikhale et al. 2020a,b,c; Patil et al. 2020; Khanal et al. 2020a,b, 2021a,b). The earlier reports indicated the Janus kinase/signal transducer and activator of transcription (JAK-STAT) pathway's deregulation during host invasion of coronavirus by affecting the type I IFNs synthesis (Luo et al 2020). But, the stimulation of IFN by host cells helps in immune modulation and prevent viruses' spreading (Li et al. 2020). Furthermore, comorbid patients with compromised immunity are more prone and carry a higher risk for COVID-19 infection (Khanal et al 2020a,b). The manipulation of the JAK-STAT pathway due to its close association with the immune system can manage the COVID-19 infection (Villarino et al. 2017; Satarker et al. 2020). Hence, the JAK-STAT pathway's multiple protein expression modulations would regulate IFN-associated genes' production (Satarker et al. 2020) to minimize viral infection.

Furthermore, coronavirus is composed of multiple proteins such as 3C-like protease (3CLpro), papain-like protease (PLpro), and spike protein that play an essential role in host entry and replication (Bhoj and Chen 2009; Lindner et al. 2005; Patil et al. 2020). Furthermore, 3CLpro alters the ubiquitin system and also modifies the functional protein by incorporating the viral polypeptides (Bhoj and Chen 2009), whereas PLpro is reported to process pp1a and pp1ab into the replicase proteins and plays a crucial task in the viral life cycle (Lindner et al. 2005). Likewise, spike protein utilizes angiotensin-converting enzyme 2 as a receptor to enter inside the cell (Kuhn et al. 2004; Patil et al. 2020). Hence, these three proteins of coronavirus could be the targets of interest for investigations to manage COVID-19.

Indian traditional system of medicines described many anti-viral plants that are effective in controlling viral infection. Given the scarcity of treatment to combat COVID-19, the present study focused on identifying JAK-STAT pathway regulators by network pharmacology and predicting the probable lead hits to act over 3CLpro, PLpro, and spike protein using in-silico molecular docking. The complete flow line for the present study is depicted in Fig. 1.

Fig. 1
figure 1

Flowline of the study to identify the JAK-STAT modulators and 3CLpro, PLpro, and spike protein inhibitors from Indian medicinal plants

Full size image

Methodology

Compounds mining and their gene expression study

The traditional anti-viral plants were identified from the open-source database and published literature. The list of each phytoconstituent was retrieved using the Chemical Entities of Biological Interest (ChEBI) database (https://www.ebi.ac.uk/chebi). Their simplified molecular input line entry system (SMILES) was queried in DIGEP-Pred (Lagunin et al. 2013) to identify the probable expression of genes based on mRNA expression at the pharmacological activity (Pa) > 0.5. A complete datasheet comprising the list of plants, their phytoconstituents, and their expression profile was constructed.

Identification of proteins involved in the JAK-STAT pathway

We identified the proteins involved in the JAK-STAT signaling pathway from the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database (pathway ID: hsa04630) (https://www.genome.jp/kegg/pathway.html). Furthermore, “Pivot Table” was used to determine the phytoconstituents and the highest count of proteins involved in the JAK-STAT pathway. Multiple keywords “STAM, JAK, TC-PTP, SHP1, STAT, SHP2\GRB, PI3K, AKT, SOS, IRF9, Ras, Raf, mTOR, PIAS, SLIM, CBPP300, CIS, Bcl-2, Bcl-XL, c-Myc, p21, AOX, GFAP, SOCS, MCL1, PIM1, and CycD” were used to identify the modulators of the JAK-STAT pathway concerning KEGG pathway; hsa04630.

Enrichment analysis, network construction, and hit identification

We performed an enrichment analysis of the JAK-STAT pathway's expressed targets by querying a list of modulated proteins using a search tool for the retrieval of interacting genes/proteins (Szklarczyk et al. 2019) and identified biological processes, cellular components, and regulated pathways concerning KEGG. The network interaction of phytoconstituents and modulated proteins from JAK-STAT pathway was constructed using Cytoscape (Shannon et al. 2003) version 3.7.1. The constructed network was analyzed based on edge count by treating the network as directed and setting the node size as “low values to small sizes” and color as “low values to bright colors”. Then the protein targeted by the highest number of phytoconstituents was further chosen as a target for docking.

Druglikeness and absorption, distribution, metabolism, excretion, and toxicity (ADMET) profile of lead hits

We predicted the probable oral bioavailability and a drug-likeness score of hit phytoconstituents using Lipinski's rule of five (Lipinski 2004). According to Lipinski's rule of five, only two violations may be accepted for human intestinal absorptivity. Similarly, the ADMET profile of lead hits was calculated using admetSAR2.0 (Yang et al. 2019). Using admetSAR2.0, multiple parameters for ADMET profile like human intestinal absorption, Caco-2, and blood-brain barrier permeability, inhibition/substrate of P-glycoprotein, CYP3A4, CYP2C9, CYP2D6, CYP3A4, CYP2C19, CYP1A2, and toxicity effects like eye corrosion/irritation, ames mutagenesis, human either-a-go-go inhibition, micronuclear, hepatotoxicity, fish aquatic toxicity, acute oral toxicity, and Tetrahymena pyriformis and binding affinity for multiple proteins like estrogen, androgen, thyroid, glucocorticoid, and aromatase receptor, PPAR-γ and plasma protein binding was predicted.

In silico molecular docking

The crystallographic protein of Tumor necrosis factor receptor-associated factor 5 (TRAF5) for humans with complete amino acid residues was not found in Research Collaboratory for Structural Bioinformatics (RCSB) protein databank (https://www.rcsb.org/). Hence, homology modeling was performed for TRAF5 using the FASTA sequence of AAH29600.1 as a query sequence against TRAF5 from Mus musculus as a template (PDB: 4GJH) using SWISS-MODEL (Schwede et al. 2003). Pre-complexed hetero-molecules from the retrieved target were removed using Discovery studio (Dassault Systèmes BIOVIA) 2019 and saved in.pdb format. The corresponding 3D structure of each ligand against TRAF5 was retrieved from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/) in.sdf format, converted into.pdb format using Discovery studio 2019; drawn in Marvin sketch (https://chemaxon.com/products/marvin) and converted into 3D.pdb format if not available. Each molecule's energy was minimized using mmff94 forcefield (Halgren 1996) and converted into.pdbqt format. Autodock4 (Morris et al. 2009) was used for docking the identified ligands with TRAF5 to get ten different ligand poses. The minimum binding energy scored pose was chosen to visualize the ligand–protein interaction using Discovery studio. Furthermore, we assessed the binding affinity of lead hits with 3CLpro (PDB: 6LU7), PLpro (PDB: 4M0W), and spike proteins (homology modeled target, accession number: AVP78042.1 as query sequence and PDB: 6VSB as a template using SWISS-MODEL) as previously explained.

Phylogeny comparison of 3CLpro, PLpro, and spike protein

The FASTA sequence of 3CLpro (accession:1P9U_A, GI: 31,616,013), PLpro (accession: 4WUR_A, GI: 726,969,143), and spike protein (accession: AVP78042.1 GI: 1,369,125,431) were retrieved from NCBI-protein and compared with the RCSB protein bank database to identify the possible homologous proteins using BlastP (Altschul et al. 1997).

Results and discussion

Identification of hits  as the JAK-STAT regulators

Twenty-five different medicinal plants were screened and seventeen possessing anti-viral activity (Abutilon indicum, Crescentia alata, Caesalpinia bonduc, Datura metel, Durio zibethinus, Evolvulus alsinoides, Garcinia mangostana, Indigofera tinctoria, Jatropha curcus, Leucas aspera, Orobanche corymbosa, Ricinus communis, Trichodesma indicum, Vitex trifolia, Wedelia Chinensis, Woodfordia fruticosa, and Wrightia tinctoria) were chosen for further study. A total of 174 multiple phytoconstituents were retrieved from the ChEBI database. Among them, 96 bioactives were identified to modulate the proteins involved in the JAK-STAT pathway. Likewise, 174 molecules were identified to regulate 807 genes; among them, 15 were from the JAK-STAT signaling pathway. Similarly, network analysis identified β-maaliene as a hit molecule to regulate the highest number of genes, and TRAF5 was a majorly targeted protein (Fig. 2). Similarly, gene ontology analysis of modulated proteins in the JAK-STAT pathway (Fig. 3) reflects the highest gene regulation in the cellular components was from cytoplasm and regulation of biological process and response to stimulus were the major biological processes.

Fig. 2
figure 2

Interaction of phytoconstituents with proteins involved in the JAK-STAT pathway

Full size image
Fig. 3
figure 3

Gene ontology analysis of regulated proteins from the JAK-STAT pathway

Full size image

Enrichment analysis of JAK-STAT modulated proteins

Enrichment analysis of proteins identified multiple pathways such as JAK-STAT signaling pathway, viral carcinogenesis, chronic myeloid leukemia, ErbB signaling pathway, apoptosis, Hepatitis B, Necroptosis, HTLV-I infection, NF-kappa B signaling pathway, Th17 cell differentiation, autophagy—animal, mTOR signaling pathway, chemokine signaling pathway, and Ras signaling pathway involved in the regulation of viral infections (Table 1).

Table 1 Enrichment analysis of modulated proteins from JAK-STAT pathway
Full size table

Druglikeness and ADMET profile of lead hits

Vitexilactone, among 39 different phytoconstituents targeting TRAF5, was predicted for the highest druglikeness score, i.e., 0.88. Among them, 18 compounds possessed a non-druglikeness character based on the rule of five. The druglikeness score of each phytoconstituent targeting TRAF5 is summarized in Table 2. Similarly, the ADMET profile of each phytoconstituent is depicted in Fig. 4 as a heat map.

Table 2 Druglikeness of phytoconstituents targeting TNF receptor-associated factor 5
Full size table
Fig. 4
figure 4

Absorption, distribution, metabolism, excretion, and toxicity profile of phytoconstituents. a (S)-3,7-dihydroxychroman-4-one, b 1-O-acetyl-4R,6S-britannilactone, c caesalpinin F, d isoceanothic acid, e jerantinine F, f mangiferin, g norcaesalpinin A, h Norcaesalpinin B, i norcaesalpinin D, j norcaesalpinin E, k platanic acid, l rotundifuran, m salicifoliol, n stigmastane-3β,5α,6β-triol, o taxine B, p vicenin-3, q viteagnusin I, r vitetrifolin D, s Vitex norditerpenoid 2, and t vitexilactone

Full size image

In silico molecular docking

Sesaminol 2-O-β-D-gentiobioside, among 39 different phytoconstituents targeting TRAF5, was predicted to have the highest binding affinity (binding energy  – 8.6 kcal/mol) via 5 hydrogen bond interactions with Gln431, Phe429, and Trp408 amino acids. However, vicenin-3 was predicted to possess the highest number of hydrogen bond interactions, i.e., 9 with Lys500, Asp502, Ser505, Ser506, Ser507, Gly520, and Ser519 amino acids; though scored comparatively lower binding affinity (binding energy  – 6.6 kcal/mol) than sesaminol 2-O-β-D-gentiobioside (Table 3). The interaction of sesaminol 2-O-β-D-gentiobioside and vicenin-3 is represented in Fig. 5.

Table 3 Binding energy and mode of interactions of phytoconstituents with TNF receptor-associated factor 5
Full size table
Fig. 5
figure 5

Interaction of a sesaminol 2-O-β-D-gentiobioside and b vicenin-3 with TNF receptor-associated factor 5

Full size image

Likewise, sesaminol 2-O-β-D-glucoside was predicted to possess the highest binding affinity (binding energy  – 9.2 kcal/mol) with 3CLpro via 2 hydrogen bond interactions with Glu166, and Thr190. Similarly, vicenin-3 was predicted to have a binding affinity (binding energy -8.2 kcal/mol) with 3CLpro via 5 hydrogen interactions with Gly143, Ser144, and Glu166 (Table 4). The interaction of sesaminol 2-O-β-D-glucoside and vicenin-3 is represented in Fig. 6.

Table 4 Binding energy and mode of interaction of phytoconstituents with 3CLpro
Full size table
Fig. 6
figure 6

Interaction of a sesaminol 2-O-β-D-glucoside and b vicenin-3 with 3CLpro

Full size image

Sesaminol 2-O-β-D-gentiobioside was predicted to possess binding affinity (binding energy  – 8.5 kcal/mol) with PLpro via 4 hydrogen bond interactions with Thr75 and Asp77. However, taxine B was predicted to have the highest hydrogen bond interactions, i.e., 6 with Thr171, Arg167, Glu204, Met207, and Gln233 amino acids though it scored comparatively lower binding affinity (binding energy -6.9 kcal/mol) with PLpro (Table 5); the interaction is represented in Fig. 7.

Table 5 Binding energy and mode of interaction of phytoconstituents with PLpro
Full size table
Fig. 7
figure 7

Interaction of a sesaminol 2-O-β-D-gentiobioside and b taxine B with PLpro

Full size image

In addition, sesaminol 2-O-β-D-gentiobioside was predicted to have the highest binding affinity (binding energy  – 9.7 kcal/mol) with spike protein via 8 hydrogen bond interactions with Thr47, Leu48, Gln825, Ala801, Ala803, Phe805, Lys807, and Asp811 (Table 6); the interaction is represented in Fig. 8.

Table 6 Binding energy and mode of interaction of phytoconstituents with spike protein
Full size table
Fig. 8
figure 8

Interaction of sesaminol 2-O-β-D-gentiobioside with spike protein

Full size image

Phylogeny comparison of 3CLpro, PLpro, and spike protein

The protein blast identified 3 different proteins, i.e., 3C-like proteinase of porcine transmissible gastroenteritis coronavirus strain Purdue, replicase of transmissible gastroenteritis virus, and Mpro Transmissible gastroenteritis virus; based on query cover (100%), E-value (0), and percentage identity (100%), these proteins were found to be similar with 3CLpro. The alignments were also matched with other proteins like Pedv 3c-like protease of diarrhea virus, replicative polyprotein 1ab of murine hepatitis virus strain A59, and purine nucleoside phosphorylase of Toxoplasma gondii. However, maximum matching was with multiple strains of the coronavirus from different organisms. Similarly, the sequence of PLpro was similar to previously reported papain-like protease from coronavirus in England and Jordan. Furthermore, PLpro sequence also matched with papain-like protease from other viral strains (Murine hepatitis virus strain A59) by 29.45%, replicase polyprotein 1ab (Murine hepatitis virus strain A59) by 29.13%, replicase polyprotein 1ab (Avian infectious bronchitis virus; strain Beaudette) by 27.11%, and nonstructural protein 3 (Avian infectious bronchitis virus; strain Beaudette) by 26.76%. Likewise, spike proteins from other viruses like Murine hepatitis virus strain A59, feline infectious peritonitis virus, infectious bronchitis virus, and porcine epidemic diarrhea virus CV777 were similar to coronavirus. Similarly, 47.06% of spike protein amino acid sequence matched with Hepcidin; however, this match's E-value was 5. The dendrogram (Fig. 9) represents the hierarchical clustering relationship of 3CLpro, PLpro, and spike protein with respective proteins.

Fig. 9
figure 9

Hierarchical clustering relationship of a 3CLpro, b PLpro, and c spike protein with multiple proteins

Full size image

The present study attempted to identify the potential lead hits of JAK-STAT signaling pathway regulators from anti-viral Indian traditional medicinal plants and identify the probable homologous proteins of COVID-19 3CLpro, PLpro, and spike protein concerning an available database. Furthermore, the study also predicted the binding affinity of a lead hit as JAK-STAT regulator with 3CLpro, PLpro, and spike protein using in silico molecular docking.

COVID-19 is reported to impair pulmonary gas exchange, primarily in the malfunctioning of the maladjusted immune system. Furthermore, to eliminate the coronavirus infection, the subject's immune system plays an important role (Li et al. 2020). KEGG database also records the JAK-STAT pathway for multiple immunological disorders, including JAK-STAT signaling for viruses (pathway entry: hsa04630). Hence, identifying lead hits to modulate the JAK-STAT pathway proteins from the medicinal plants can play an essential role in regulating an immune system as a home remedy and may help as a preventive approach to deal with corona infection. Furthermore, we were interested in knowing if 3CLpro, PLpro, and spike protein share their amino acid sequence with any other predeposited proteins in the RCSB databank.

Interestingly, we identified 3CLpro to share its amino acid sequences with other microbial infections like PEDV main protease from Porcine epidemic diarrhea virus CV777. Furthermore, PLpro also shared its amino acid sequence with replicase polyprotein 1ab of Murine hepatitis virus strain A59. Likewise, spike protein possessed similar amino acid alignments with spike glycoprotein of Murine hepatitis virus strain A59. This observation suggests that analogues of protease inhibitors from procaine epidemic diarrhea could be chosen as the investigative molecule against 3CLpro. Similarly, analogues of replicase polyprotein 1ab and spike glycoprotein of Murine hepatitis virus strain A59 could act over the PLpro and spike protein, respectively.

Enrichment analysis of modulated proteins from JAK-STAT pathways identified the multiple pathways which may contribute to maintain an immune system like apoptosis, necroptosis, Th1, and Th2 cell differentiation, NF-kappa B signaling pathway, Th17 cell differentiation, and chemokine signaling pathway. Coronavirus replication is linked with the membrane of the endoplasmic reticulum and is reported to increase stress over it; leads to apoptosis via unfolded protein response (Fung and Liu 2014). Likewise, the corona virus-induced necrotic cell death may occur due to multiple mechanisms like Rip3-induced oligomerization, mixed lineage kinase domain-like protein dispensed cell death, necroptotic elements-induced human lung cell damage, lysosome-induced damage, and activation of multiple inflammasomes (Yue 2018). Hence, modulation of these pathways’ may contribute to the immune system’s homeostatic regulation and reduce cell apoptosis and necrosis. Other pathways involving T-helper cells, cytokine, and chemokine signallers are also closely associated in enhancing the immune system that may contribute towards up-/down- regulation of the immune system during corona infection (Padovan 2017; Simon 2011; Jang et al. 2005).

The network pharmacology screening of multiple phytoconstituents from folk medicines has identified the β-maaliene to regulate the highest number of proteins. The probable checkpoints affected by β-maaliene in the JAK-STAT pathway are presented in Fig. 10. Although, based on the literature search, we did not find any reported anti-viral activity of β-maaliene, our prediction of the biological spectrum via Prediction of Activity Spectra for Substances (PASS) ONLINE (Lagunin et al. 2000) identified it to possess an anti-viral effect against influenza, adenovirus, Cytomegalovirus, Herpes, Picornavirus, and Hepatitis C. Furthermore, the compound was also envisioned as viral entry inhibitor. So, apart from modulating the JAK-STAT pathway proteins, β-maaliene may also act over the spike protein of coronavirus and inhibit its entry into the cell. Furthermore, these calculations suggest that analogs of β-maaliene could be the choice of investigative molecules as a modulator of the JAK-STAT pathway and inhibitor of virus entry into host cells.

Fig. 10
figure 10

Probable checkpoints affected by β-maaliene in JAK-STAT pathway. ★targeted proteins by β-maaliene

Full size image

We identified TRAF5 protein that may be modulated by the highest number of phytoconstituents during network pharmacology-evident target identification. A previous report suggested that the requirement of TRAF5 optimizes the T cell expansion in response to infection (Kraus et al. 2008). Hence, we attempted to identify the lead hit agent as the prime modulator of TRAF5 from the list of phytoconstituents. Based on the binding energy and number of hydrogen bond interactions we identified sesaminol 2-O-β-D-gentiobioside and vicenin-3 to possess the best interaction with TRAF5. Furthermore, TRAF5 modulators were also docked with 3CLpro, PLpro, and spike protein; and found that sesaminol 2-O-β-D-gentiobioside, taxine B, vicenin-3, and sesaminol 2-O-β-D-glucoside as lead hits. In addition, PASS (Lagunin et al. 2000) helps to identify the probable anti-viral compounds by querying the SMILES of molecules which predicted plausible anti-viral properties against Influenza (Pa: 0.733), Hepatitis B (Pa: 0.402), and Herpes (Pa: 0.371) by sesaminol 2-O-β-D-gentiobioside (1), Rhinovirus (Pa:0.324) by Taxine B (2), Herpesvirus (Pa: 0.842), Influenza (Pa: 0.540), Hepatitis B virus (Pa: 0.443), Picornavirus (Pa: 0.349), Poxvirus (Pa: 0.237), Hepatitis (Pa: 0.136), and Retrovirus (Pa: 0.133) by vicenin-3 (3) and Influenza (Pa: 0.690), Hepatitis B (Pa: 0.378), Herpes (Pa: 0.396), Rhinovirus (Pa:0.305) by sesaminol 2-O-β-D-glucoside (4). Docking study also revealed the probable binding affinity of these molecules with 3CLpro, PLpro, and spike protein. These targets also shared the similar amino acid sequence with proteins from other viruses; hence, these molecules may be further screened via their functional group modeling to target coronavirus with the proper assessment of protein domain of target of interest.

Although the present screening focused in identifying the lead hit anti-viral molecules against COVID-19, Indian herbal medicinal plants were screened via network pharmacology and gene enrichment analysis; is a concept via "multi compounds-multi protein" interaction. Furthermore, the present study's outcome is knowledge-based and depends on multiple regression models for particular biological activities. On the other hand, during COVID-19 infection, the dysregulated immune system may play an essential role in spreading the virus in the host system; it could be maintained via the JAK-STAT signaling pathway, the main basic concept of the present study. Likewise, identifying the homologous proteins may also play a crucial role in identifying new drug molecules that kindled us for phylogeny comparison of targets based on the predeposited amino acid sequences of individual proteins.

Conclusion

The present study identified few lead hits from the traditional Indian medicines reported for their anti-viral property against COVID-19. Further phylogeny comparison also reflected that the coronavirus proteins may share a common domain with other similar microbial proteins. It also suggests that reported inhibitors of such proteins could be implemented in screening the molecules as 3CLpro, PLpro, and spike protein inhibitors. However, the present findings are based on available databases and computer simulations. Hence, further investigations using wet-lab protocols are required to confirm the conclusions accordingly.