Integrative multiomics and in silico analysis revealed the role of ARHGEF1 and its screened antagonist in mild and severe COVID-19 patients.

COVID-19 is a sneaking deadly disease caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). The rapid increase in the number of infected patients worldwide enhances the exigency for medicines. However, precise therapeutic drugs are not available for COVID-19; thus, exhaustive research is critically required to unscramble the pathogenic tools and probable therapeutic targets for the development of effective therapy. This study utilizes a chemogenomics strategy, including computational tools for the identification of viral-associated differentially expressed genes (DEGs), and molecular docking of potential chemical compounds available in antiviral, anticancer, and natural product-based libraries against these DEGs. We scrutinized the messenger RNA expression profile of SARS-CoV-2 patients, publicly available on the National Center for Biotechnology Information-Gene Expression Omnibus database, stratified them into different groups based on the severity of infection, superseded by identification of overlapping mild and severe infectious (MSI)-DEGs. The profoundly expressed MSI-DEGs were then subjected to trait-linked weighted co-expression network construction and hub module detection. The hub module MSI-DEGs were then exposed to enrichment (gene ontology + pathway) and protein-protein interaction network analyses where Rho guanine nucleotide exchange factor 1 (ARHGEF1) gene conjectured in all groups and could be a probable target of therapy. Finally, we used the molecular docking and molecular dynamics method to identify inherent hits against the ARHGEF1 gene from antiviral, anticancer, and natural product-based libraries. Although the study has an identified significant association of the ARHGEF1 gene in COVID19; and probable compounds targeting it, using in silico methods, these targets need to be validated by both in vitro and in vivo methods to effectively determine their therapeutic efficacy against the devastating virus.

INTRODUCTION

Eighteen years after the outbreak of severe acute respiratory syndrome (SARS) caused by SARS‐coronavirus (SARS‐CoV), the world is now facing a new dare posted by version 2.0 of CoV causing similar symptoms like SARS emerged from Wuhan city, China, which was initially named SARS‐CoV‐2 and later renamed by the World Health Organization as COVID‐19 on February 11, 2020. , ,  As of July 1, 2021, the total number of laboratory‐confirmed cases of COVID‐19 infection in the world were 182 988 210, including 3 962 976 associated deaths. The first COVID‐19 patient‐reported from Seafood Market in Wuhan was possibly infected by the direct or indirect transmission of COVID‐19 virus from bats or pangolins and later the epidemic spread to the whole world via human‐to‐human transmission. , The spread toward Asian countries, especially in India, has further raised the concern to the wider population.

CoV belongs to the largest group of viruses and comes under the Coronaviridae family and the Coronavirinae subfamily, found in birds and mammals. CoVs are classified into four genera: α, β, δ, and γ where COVID‐19 refer to β‐CoV which also consists of Middle East respiratory syndrome (MERS)‐CoV and SARS‐1. (p. 2),

Images from transmission electron microscopy revealed that CoV gives a crown‐like appearance with spike proteins projecting out from the spherical virion surface. COVID‐19 shares about 50% similarity with MERS whereas approximately 79% with SARS‐CoV indicates that SARS‐CoV‐2 is closer to SARS‐CoV. All the human CoV consists of 3ʹ‐polyadenylated tails which encodes structural proteins (nucleocapsid protein, membrane proteins, spike proteins, and envelope proteins) responsible for maintaining viral life cycle, and 5ʹ‐methylated cap region is for nonstructural proteins, which are important for replication of the virus. , (p. 19), ,

Among various human CoVs, SARS‐CoV‐2 caused a severe mortality rate close to 1.38, and the most common symptoms in positive cases were cough, fever, headache, fatigue, diarrhea, and blood‐stained mucous. In the current scenario, the characterization of novel biomarkers which may play an important role in prediction and observation in the status of SARS‐CoV‐2 disease is required. , (p. 2)

Many studies related to SARS‐CoV‐2 have been published in these 2 years covering aspects like high‐throughput transcriptomic analysis, network‐based dynamics, artificial intelligence‐based drug repurposing, and novel therapeutic compounds identification followed by molecular dynamics (MD) validation. Our study is novel as it presents an age‐based weighted network of SARS‐CoV‐2 transcriptomic dataset for robust identification of highly correlated hub module along with therapeutic targets identification thus enabling us to understand the disease mechanism at the molecular and systems level.

In the present study, SARS‐CoV‐2 associated messenger RNA (mRNA) expression profile was retrieved from National Center for Biotechnology Information (NCBI)–Gene Expression Omnibus (GEO) followed by identification of differentially expressed genes (DEGs) in different infection severity groups. The overlapping mild and severe infectious (MSI)‐DEGs present in both groups were subjected to trait‐linked weighted gene co‐expression network (WGCN) construction and hub module detection. The hub module MSI‐DEGs were forwarded to protein–protein interaction (PPI) network construction and enrichment analyses for predicting potential therapeutic targets. Lastly, we applied molecular docking and MD method against Rho guanine nucleotide exchange factor 1 (ARHGEF1) gene to discover potential hits for SARS‐CoV‐2 from antiviral, anticancer, and natural product‐based libraries. The novelty of our work involves an amalgamation of transcriptomic, network analysis, PPI, and enrichment analyses of COVID‐19 peripheral blood mononuclear cell (PBMC) patient samples for investigation of human target gene which might be responsible for controlling the severity of SARS‐CoV‐2 infection, unlike earlier research work. , The potential candidate hits can be further selected for in vitro and in vivo studies in future work to speed up the drug discovery against SARS‐CoV‐2.

MATERIALS AND METHODS

SARS‐CoV‐2 microarray data collection and differential expression analysis

NCBI‐GEO (https://www.ncbi.nlm.nih.gov/geo/) was queried to extract SARS‐CoV‐2‐related mRNA expression profile with “COVID‐19” and “SARS‐CoV‐2” being used as appropriate keywords. Our searching criteria for picking a suitable dataset was as follows: (1) the dataset must be “expression profiling by array” type and its samples must be from “Homo Sapiens”; (2) the dataset must be having both healthy and COVID‐19 (varying severity levels) patient samples; (3) the dataset must be having diverse age group samples; (4) the dataset must be having processed and raw microarray data; and (5) the dataset submission date must be within 1 year. Any studies devoid of nonhuman samples, case reports, review articles, abstracts, and cell‐line‐based experimental study designs were excluded. Series matrix expression file of the chosen dataset was downloaded using GEO query R package. Mapping of probe IDs to their corresponding HUGO Gene Nomenclature Committee (HGNC) symbols was performed using feature data stored in the expression set object of dataset. The averaging expression value of genes mapping to more than one probe IDs was done to avoid redundancy. The COVID‐19 patient samples in the dataset were bifurcated into groups based on their infection severity levels followed by the elimination of 50% low variance genes. HGNC multisymbol checker (https://www.genenames.org/tools/multi-symbol-checker/) was utilized and only officially approved HGNC symbols were retained. A two‐sample statistical t‐test was utilized for calculating the p‐value of each gene between control and COVID‐19 samples (infection severity group‐based) followed by obtaining their and Benjamini–Hochberg (BH) p‐values utilizing the Limma R package. The genes corresponding to BH p < 0.01 and were regarded as DEGs. DEGs with and were bifurcated as up‐ and downregulated, respectively. The DEGs overlapping between different infection severity groups was considered as MSI‐DEGs.

WGCN construction and hub module identification

Weighted Gene Co‐Expression Network Analysis (WGCNA) R package was utilized for establishing a WGCN from SARS‐CoV‐2 specific MSI‐DEGs and identifying representative module genes having a high correlation with clinical characteristics (i.e., age). All the MSI‐DEGs along with their samples were passed through the goodSamplesGenes function to eliminate any missing values. The samples were clustered thereafter to eliminate any outliers. The age information of all samples was also taken into consideration before identifying modules. The pickSoftThreshold function assisted in selecting a suitable soft‐thresholding power () to which co‐expression similarity will be raised for computing adjacency. In general, a suitable value of was chosen in compliance with the approximate scale‐free topology criterion. However, in some cases, the data does not follow appropriate scale‐free topology due to the presence of a strong driver which makes a subgroup of samples globally distinct from the rest. This difference causes an elevated correlation amongst a large cluster of genes that invalidates the assumption of the scale‐free topology approximation. In such cases, a suitable value of needs to be chosen based on the sample size in accordance with the WGCNA guidelines (https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/faq.html) to make the resulting network conservative. For reduction in noise and false associations, the weighted adjacency matrix was converted into a topological overlap matrix (TOM) followed by computation of corresponding dissimilarity (dissTOM). The hclust function was utilized to generate a hierarchical clustering tree (dendrogram) of genes in consideration with the dissTOM measure. To identify densely interconnected highly co‐expressed gene patterns (i.e., modules) from the branches of a tree, a dynamic tree cut algorithm was incorporated. Module eigengene (ME) and dissimilarity measure between MEs were computed to merge the modules with highly co‐expressed genes. ME dendrogram was checked based on Pearson correlation for merging multiple modules with comparable expression profiles. Correlation‐based absolute module significance (i.e., average gene significance [GS] of participating genes in a given module) with the age of samples followed by module membership (MM) (correlation of the ME and the gene expression profile) for all modules were computed. Intramodular connectivity (k.in) can be explained as a measure of network connectivity with respect to nodes or genes of a specific module. The correlation of MM versus GS, GS versus k.in, and MM versus k.in was used to identify the most significant associations and the module having significantly highest correlation with age was chosen as the hub module.

PPI network construction, pathway, and gene ontology term enrichment analyses

The MSI‐DEGs present in our hub module were subjected to PPI network construction using Search Tool for the Retrieval of Interacting Genes (STRING, https://string-db.org/) v11.5 web‐based tool. This PPI network was formed at the highest confidence (corresponding to interaction score ) and afterward visualized using Cytoscape v3.9.0. Pathway and GO term enrichment data for MSI‐DEGs present in the hub module was compiled using various libraries (i.e., reactome for pathways and gene ontology [GO]‐biological process [BP], GO‐molecular function [MF], and GO‐cellular compartment terms) available within Enrichr database. Top 10 pathways and GO terms corresponding to p  were regarded as statistically significant ones. The MSI‐DEGs overlapping between constructed highest confidence PPI network, top 10 significant pathways, and GO terms were regarded as SARS‐CoV‐2 hub MSI‐DEG(s), respectively.

Identification of binding sites and preparation of protein

RhoGEF1 enzyme encoded by ARHGEF1 gene was the hub driver MSI‐DEG in SARS‐CoV‐2 infectious process and was selected as a promising target for COVID‐19. The binding sites of the identified protein (RhoGEF1), PDB ID: 3ODO, were searched using the machine learning method on PrankWeb server (https://prankweb.cz/) due to very little information available about protein binding site. Before this, the protein structure was prepared using the automatic protein preparation module of Discovery Studio (DS) 2019. Hydrogen atoms were added followed by side‐chain refinement and protonation states were determined. A receptor grid sphere with coordinates X = 46.5266, Y = −41.5728, Z = 26.6118, and radii = 12 Å was generated based on PrankWeb server information.

High‐throughput virtual screening using LibDock

LibDock is a high‐throughput screening program in which both ligands, as well as protein, were rigid. It enumerates the hotspots at the binding site of the protein by using apolar and polar probes which were used to match the ligands to form favorable interactions. Natural‐Product‐Like Compound Library (900 compounds), Anticancer Screening Library (6300 compounds), and Antiviral Compound Library (13 700 compounds) were retrieved from Life Chemicals (https://lifechemicals.com/) and COVID‐19 related library (35 compounds) to screen ARHGEF1 inhibitors. Using the defined coordinates site, an active site for molecular docking was generated. Structure‐based virtual screening was executed by docking all the compound libraries at the prepared active site by using the LibDock module of DS 2019. The LibDock score was used as a criterion to rank the docked poses and grouped by name.

Ligand preparation and ADMET prediction

The filtered ligands from LibDock docking were selected for the calculation of the ADMET properties (absorption, distribution, metabolism, excretion, and toxicity) and Ames mutagenicity using ADMET and TOKAT module of DS 2019. These compounds were further filtered based on their molecular properties using Lipinski's and Veber rule implemented in DS 2019. Filtered compounds were subjected to energy minimization by using the steepest descent (5000 cycles) followed by the conjugate gradient (5000 cycles) to satisfy the minimum criteria of root mean square gradient of 0.001.

Molecular docking

All the ligands that were successfully filtered out were submitted for molecular docking using the CDOCKER module of DS 2019 which is a CHARMM‐based MD simulated docking tool in which ligands were allowed to be flexible while receptor to remain rigid. Crystal structure of ARHGEF1 (PDB ID: 3ODO) with the same binding site was selected for molecular docking. By retaining all the default parameters, filtered ligands were docked into the binding site of ARHGEF1. Different conformations for each ligand were generated and investigated based on CDOCKER energy. Next, we generated the nonbond interactions between the ARHGEF1, and ligand conformations obtained from CDOCKER docking with the help of Analyse ligand poses module of DS 2019. Detailed knowledge of nonbond interactions is important in structure‐based design, as it helps to select the favorable interactions between ligands and their receptor.

Flexible docking

To cross‐validate the identified binding site and molecular docking process, a flexible docking method was used. It may also generate more correct or accurate ARHGEF1 complexes which may be used as a starting structure for MD simulation. The top 15 complexes were selected after analyzing the nonbond interactions for redocking by using the flexible docking module of DS 2019 which docks the hit ligands with induced‐fit protein optimization while simulating the protein flexibility. Favorable residues from analyzing ligand pose protocol were defined as flexible residues. The receptor conformations (50) were generated using ChiFlex and subsequently, 255 conformations were generated for ligand. Next, ligand docking was done by the LibDock algorithm followed by side‐chain refinement using ChiRotor. The final refinement of complexes was done by using CDOCKER under the CHARMM forcefield.

MD study

The top three hit‐complexes were used as starting structures for MD simulation by using Gromacs 4.5.6 package and Gromos96 forcefield. Additionally, the top two hit complexes from the COVID‐19 related library were also included in this study. The topology and parameter file for ligands were generated on the GlycoBioChem PRODRG2 Server (http://davapc1.bioch.dundee.ac.uk/cgi-bin/prodrg). The complexes were solvated with a cubic box of a simple point charge with 1 Å margin. Systems were next neutralized by adding Na+ and Cl− ions to balance the system charges. Each complex was energy minimized with 10 000 steps of the steepest descent algorithm to remove any steric clashes. Then the system was equilibrated with 100 ps of NVT followed by 100 ps of NPT at a constant temperature of 300 K and a constant pressure of 1 atm. All bonds were constrained with LINC algorithm and long‐range electrostatics were calculated with the particle‐mesh Ewald algorithm. Finally, 50 ns of MD run was performed and snapshots were saved for every 5 ps. Final trajectories were analyzed in pymol and Gromacs tools (g_energy, g_rmsf, g_rmsd, g_gyrate, and g_hhbond) to calculate potential energy, root‐mean‐square fluctuation (RMSF), root‐mean‐square deviation (RMSD), the radius of gyration, and h‐bond analysis.

Molecular mechanics Poisson–Boltzmann surface area

The binding free energies for docked complexes were calculated using molecular mechanics Poisson–Boltzmann surface area (MM‐PBSA) to revalidate the MD simulation studies using the g_mmpbsa tool. We performed the calculations from the last 5 ns of each MD run and used all the set parameters. Next, the binding energies of the best three ligands‐complexes were enumerated using the following equations as reported earlier :

RESULTS

On the basis of searching criteria, we chose SARS‐CoV‐2 mRNA expression profile possessing accession number GSE164805. It comprised a total of 15 PBMC patient samples (i.e., 5 controls and 10 COVID‐19 [5 mild + 5 severe]). Mapping probe IDs to their corresponding genes was followed by averaging expression values corresponding to duplicates which led to 44 822 unique genes. The COVID‐19 patient samples in the dataset were bifurcated into MSI groups along with common controls in both groups leading to 10 samples in each group (i.e., Group 1: 5 mild + 5 controls and Group 2: 5 severe + 5 controls). The variance of each gene across samples in both groups was computed leaving 22 411 high variance genes in both the groups. After passing these genes through HGNC multisymbol checker, a total of 10 729 and 10 518 officially approved HGNC symbols were left in Groups 1 and 2, respectively. A total of 837 and 1935 genes were differentially expressed corresponding to BH p  and in Groups 1 and 2. Amongst these DEGs, 427 and 1129 were downregulated while 410 and 806 were upregulated in Groups 1 and 2. Volcano plots showing the distribution of significant (colored points) as well as nonsignificant (black colored points) genes in Groups 1 and 2 were shown in Figure 1A,C. Figure 1B,D shows an annotation heatmap of top 10 up and downregulated DEGs in Groups 1 and 2. The sample type, gender, and age annotation bars were placed at the top of the heatmap. The controls clustered distinctly from mild and severe samples in both groups as evidenced from the heatmap. FST [] and TEX101 [] were most down and upregulated in Group 1. Whereas ADAT1 [] and TEX101 [] were most down and upregulated in Group 2. Among the top 10 up‐ and downregulated DEGs, ADAT1, HBD, and FST were commonly downregulated in both groups. Also, TEX101, FAM177B, and OR51E1 were commonly upregulated in both groups. TEX101 was the highest upregulated DEG in both groups. Group 1 samples constituted 20% females and 80% males while Group 2 samples constituted 10% females and 90% males. Venn plot showing the 569 overlapping MSI‐DEGs between both groups is shown in Figure S1A. Figure S1B shows a kernel density plot distribution of sample age for all 569 MSI‐DEGs. Three samples (i.e., one control, one mild, and one severe) with age 54 were having the highest density as evidenced from the plot. All the samples in Group 2 were having an age with 73 being the maximum age. Whereas all samples in Group 1 except one were had age with 71 being the maximum age.

Figure 1

(A) Volcano plot showing the distribution of DEGs (downregulated: 427 +  upregulated: 410) and nonsignificant (9892) genes across control and mild COVID‐19 samples (i.e., Group 1). (B) Annotation heatmap showing the distribution top 10 up‐ and downregulated DEGs across Group 1 samples. The expression value of DEGs (in rows) is normalized across all the samples (in columns). Cluster dendrograms representing Euclidean distance‐based hierarchical clustering for both rows and columns are presented along the left and top sides of the plot. Sample type (medium blue for controls and carmine for mild), gender (eastern blue for females and midnight blue for males), and age annotation bars are placed at the top of the heatmap. (C) Volcano plot showing the distribution of DEGs (downregulated: 1129 + upregulated: 806) and nonsignificant (8583) genes across control and severe COVID‐19 samples (i.e., Group 2). (D) Annotation heatmap showing the distribution top 10 up‐ and downregulated DEGs across Group 2 samples. The expression value of DEGs (in rows) is normalized across all the samples (in columns). Cluster dendrograms representing Euclidean distance‐based hierarchical clustering for both rows and columns are presented along the left and top sides of the plot. Sample type (Prussian blue for controls and ruby red for severe), gender (green for females and orange for males), and age annotation bars are placed at the top of the heatmap. The green‐ and magenta‐colored points signify the up‐ and downregulated DEGs in volcano plots. While the black‐colored points signify the nonsignificant genes. The x‐ and y‐axes represent (fold change) and (p‐value). DEG, differentially expressed gene

The expression data of all 569 MSI‐DEGs along with their sample age information were used as an input in WGCNA. No obvious sample outliers were detected from the sample clustering dendrogram as shown in Figure S2. Our WGCN did not follow the scale‐free topology criterion for choosing the appropriate value of with respect to . Hence, in accordance with the WGCNA guidelines, (for less than 20 samples) was chosen. Hierarchical clustering tree and dynamic tree cut algorithm revealed a total of five color‐coded gene modules (i.e., blue, turquoise, brown, green, and yellow) as shown in Figure 2A. The size of these five gene modules were as follows: blue = 109, turquoise = 231, brown = 102, green = 56, and yellow = 71. The multidimensional scaling plot of all modules in three scaling dimensions is shown in Figure 2B. There was no need for merging these modules due to the low merging height observed in the ME dendrogram. Figure S3A shows the gene co‐expression network as a heatmap plot depicting TOM among all module genes. Figure 2C shows the Barplot of GS (correlated with age) with respect to each module. The GS values for the modules were as follows: blue = 0.14, turquoise = 0.27, brown = 0.15, green = 0.41, yellow = 0.23. MM versus GS correlation values along with their p‐values for each module are shown in Table S1. The GS versus k.in values along with their p‐values for each module is shown in Table S2. MM versus k.in values along with their p‐values for each module are shown in Table S3. As evidenced from Tables S1–3, blue, yellow, and turquoise modules were eliminated from any further analysis since they were not statistically significant (). All these results indicated that the green module was the most promising and highly correlated (module significance , GS versus MM , GS versus k.in , MM versus k.in ) as compared to the brown module. Therefore, the green module (comprising 56 MSI‐DEGs) was chosen as our hub module. Figure 2D shows a scatterplot of GS for age with respect to MM in the green hub module. Figure S3B shows a scatterplot of GS for age with respect to MM in the brown module. Figure S3C,D shows heatmap plots of green and brown module genes along with their corresponding ME levels.

Figure 2

(A) Hierarchical clustering dendrogram of MSI‐DEGs clustered based on the dissimilarity measure (dissTOM) and five color‐coded modules (obtained using Dynamic Tree Cut). The modules contained highly similar expression profiles with the following sizes: turquoise (231), blue (109), brown (102), yellow (71), and green (56). (B) Three‐dimensional MDS plot where each colored point denotes a gene belonging to the module of the corresponding color. (C) Barplot showing the distribution of average gene significance (GS), that is, “module significance” and error bars across turquoise, blue, brown, yellow, and green modules. (D) Scatterplot showing significantly () highest correlation of GS for age with module membership (MM) in the green module. DEG, differentially expressed gene; MDS, multidimensional scaling; MSI, mild and severe infectious; TOM, topological overlap matrix

PPI network construction, pathway, and GO term enrichment analyses

All the 56 MSI‐DEGs present in our hub module were given as an input to the STRING web‐based tool where a total of 4 (3 upregulated + 1 downregulated) MSI‐DEGs participated in the PPI network corresponding to STRING interaction score . The PPI network as shown in Figure 3A comprised four nodes and four edges. Also, a total of 12, 14, and 10 out of 56 MSI‐DEGs within our hub module participated in the top 10 significant pathways, GO‐BP, and GO‐MF terms. The most significant pathway, GO‐BP, and GO‐MF terms were G alpha (12/13) signaling events (), regulation of mitotic spindle checkpoint (), and ATP‐dependent DNA helicase activity (). Sankey plot showing the association of top 10 significant GO‐BP and GO‐MF terms with corresponding 18 genes is shown in Figure 3B. Chord plot showing the association of 12 MSI‐DEGs with the top 10 significant pathways is shown in Figure 3C. The width of interaction edges connecting the pathways and MSI‐DEGs varied according to the p‐values. PSMF1, ARHGEF9, GNG4, and AP2M1 appeared in a maximum number of pathways (i.e., 4) whereas F10, ABCG2, BRIP1, GJD3, RASAL2, and DUSP9 appeared in a minimum number of pathways (i.e., 1). Venn plot as shown in Figure 3D revealed ARHGEF1 as the overlapping hub MSI‐DEG between PPI, pathway, and GO term (i.e., GO‐BP, GO‐MF) genesets.

Figure 3

(A) Protein–protein interaction (PPI) network comprising four nodes and four edges corresponding to STRING interaction score , that is, “highest confidence score.” The red‐ and green‐colored nodes represent up‐ and downregulated proteins. (B) Sankey plot showing the association of top 10 significant GO‐BP and GO‐MF terms with corresponding 18 genes. The width and color intensity of interaction edges varied with respect to the p‐values. (C) Chord plot showing the connection of 12 MSI‐DEGs (on right semicircle) with 10 significant pathways (on left semicircle) via colored edges. The width of edges varied with respect to p‐values. The edges initiated from unique colored strips present on the right semicircle (indicating genes) and converged to unique colored strips present on the left semicircle (indicating pathways). (D) Venn plot showing the overlapping hub MSI‐DEG (1) between PPI, GO‐BP, GO‐MF, and pathway genesets. BP, biological process; DEG, differentially expressed gene; GO, gene ontology; MF, molecular function; MSI, mild and severe infectious

Since there is no crystal structure available for our predicted target ARHGEF1 with co‐crystallized ligand and very little knowledge available about binding site, the refined or prepared structure of ARHGEF1 protein (PDB ID:3ODO) was examined on PrankWeb server and topographic features of protein were assessed through P2Rank program. Out of the three predicted pockets, we selected pocket 1 based on rank score. The amino acids predicted in pocket 1 were Thr427, Ala430, His431, Met434, Arg551, Leu552, Asp556, Met557, Thr560, Gln563, Arg564, and Lys567. On the basis of these amino acids found in the pocket (score = 5.22) (Table S4 and Figure S4), we defined the binding site for structure‐based virtual screening.

To identify novel compounds against ARHGEF1 through Dbl‐ and pleckstrin‐homology (DH/PH) binding pocket, a docking‐based virtual screening was performed using the LibDock module of DS 2019. Total 20 900 purchasable compounds comprising antiviral, anticancer, and natural products were taken from Life Chemicals (https://lifechemicals.com/) database and 35 COVID‐19 related inhibitors from Selleck chem library which were currently tested for COVID therapy (https://www.selleckchem.com/covid-19-related-products.html) were also selected as a positive control. A total of 741 542 conformations of ligand were generated which aligned to the polar or apolar region of the protein known as hotspots and 1 491 483 poses were docked successfully into the predicted binding site of the ARHGEF1. All the docked poses were shortlisted based on the highest LibDock score.

To become a successful drug, ADMET, TOPKAT, and physiochemical properties of hit‐compounds were predicted by using different modules of DS 2019. The top 300 compounds were selected for toxicity prediction followed by ADMET properties. Next, the oral bioavailability of selected hit compounds was calculated by Lipinski's rule of 5 and Veber's rule. Total 197 compounds following the TOPKAT, ADMET, Lipinski, and Veber's rule were subjected to energy minimization via steepest and conjugate algorithm. As COVID‐19 related compounds were already approved for other targets, there is no need to perform ADMET prediction. Therefore, final hit compounds were having satisfying drug‐likeness properties and were predicted to have good bioavailability.

A total of 197 filtered compounds, were docked into the predicted binding site (X = 46.5266, Y = −41.5728, Z = 26.6118, and radii = 12 Å) of ARHGEF1 protein (PDB ID:3ODO). Out of 197 compounds, 125 ligands were able to dock into the defined receptor cavity. Additionally, only 5 out of 35 COVID‐19 related compounds were able to interact with ARHGEF1 using a favorable cDocker energy score. All the ligands were arranged according to cDocker energy score whereas positive cDocker energy scored compounds were removed from the list. The molecular interaction of the binding poses was analyzed with analyse ligand poses in DS 2019. From all the binding poses, a total of 2381 favorable interactions, 1088 hydrogen bonds, 1052 hydrophobic interactions, and only 19 interactions were unfavorable. The top five residues with favorable interactions were: Lys567, Arg551, Met434, Ala430, and Leu552 (Figure S5).

Flexible docking or induced fit docking overcome the problem of false‐negative results and was more precise and accurate but it took more time to dock one ligand into the receptor cavity than other methods. The top 10 ligands from CDOCKER docking results and 5 ligands from COVID‐19 related compounds were redocked via flexible docking method into the above binding site using DS 2019. The obtained docking poses were sorted based on their LibDock score (first‐level screening) followed by cDocker energy (second‐level screening). The cDocker energy score after flexible docking (third‐level screening) of each pose was improved tremendously and enumerated in Table 1 and interactive residues of ARHGEF1 with top three hits were shown in Table 2 whereas docking scores of top hits from COVID‐19 related compounds were enumerated in Table S5 and top two ligand interactions were shown in Table S6. The receptor–ligand interactions were analyzed in DS 2019 and, for instance, the interactions of the top three ligands, F3222‐4380, F3394‐0943, and F6548‐1087 are explained in Figure 4 whereas interactions of the top two ligands (COVID‐19 related library) are explained in Figure S6. It can be seen in Figure 4A, F3222‐4380 binds with the active site of ARHGEF1 through interactions: hydrogen bonds (His431, Arg550, and Arg551), pi‐cation/alkyl/pi‐alkyl (Ala430, Arg433, Leu552, Lys567, and Arg564), and van der Waals interactions (Ser394, Thr427, Met434, Val437, Gln563, and Thr560). F3394‐0943 interacted via hydrogen bond (Thr427, Arg551, and Lys567), pi‐cation/alkyl/pi‐alkyl (Ala430, His431, Arg433, Arg550, Arg564, and Lys567), and with van der Waals interactions (Met434, Val437, Leu552, and Gln563) (Figure 4B). Finally, the interaction between F6548‐1087 and ARHGEF1 was also observed, the interactions were hydrogen bonds (His431, Arg564, and Lys567), pi‐cation/alkyl/pi‐alkyl (Ala430, Arg433, and Leu552), and van der Waals (Glu423, Thr427, Met434, Val437, Arg550, Arg551, and Thr560) (Figure 4C).

Table 1

Molecular docking score of the top 10 selected hits

Compound	Name	First‐level screening	Second‐level screening	Third‐level screening	Class
		LibDock score	‐CDOCKER energy	Flexible (‐CDOCKER energy)
Lead 1	F3222‐4380	132.739	34.7124	50.0296	Antiviral
Lead 2	F3394‐0943	135.581	30.5814	45.7904	Anticancer
Lead 3	F6548‐1087	129.865	28.3122	44.6191	Anticancer
Lead 4	F3394‐0911	135.581	30.5814	43.8092	Anticancer
Lead 5	F0372‐0536	137.649	34.5305	43.4841	Antiviral
Lead 6	F0375‐0251	138.437	26.6627	42.1999	Anticancer
Lead 7	F1342‐0042	131.347	32.7999	41.738	Antiviral
Lead 8	F3407‐3720	126.415	30.3855	41.2499	Anticancer
Lead 9	F0886‐0075	136.183	27.0587	40.27	Antiviral
Lead 10	F0886‐0050	130.1	32.6969	36.0879	Antiviral

Table 2

Nonbond interactions between RhoGef1 and top three compounds

Compound name	Residues in contact	Types of interaction	Distance (A°)
Lead 1 (F3222‐4380)	His431	Carbon hydrogen bond	2.33
	Arg550	Carbon hydrogen bond	2.76
	Arg551	Conventional hydrogen bond	3.79
	Ala430	Alkyl	4.39
	Arg433	Pi‐alkyl	4.69
	Leu552	Alkyl	4.97
	Arg564	Pi‐cation	4.58
	Lys567	Pi‐alkyl	2.55
Lead 2 (F3394‐0943)	Thr427	Carbon hydrogen bond	3.59
	Arg551	Conventional hydrogen bond	5.73
	Lys567	Conventional hydrogen bond	5.12
	Ala430	Pi‐alkyl	3.79
	His431	Pi‐alkyl	4.58
	Arg433	Pi‐alkyl	4.38
	Arg550	Pi‐alkyl	5.49
	Arg564	Pi‐cation	4.61
	Lys567	Pi‐alkyl	6.19
Lead 3 (F6548‐1087)	His431	Carbon hydrogen bond	4.77
	Arg564	Conventional hydrogen bond	4.83
	Lys567	Conventional hydrogen bond	5.54
	Ala430	Pi‐alkyl	5.02
	Arg433	Pi‐alkyl	5.25
	Leu552	Pi‐alkyl	4.97

Figure 4

2D and 3D molecular interactions between ligands and ARHGEF1obtained by molecular docking. (A) F3222‐4380, (B) F3394‐0943, and (C) F6548‐1087. All the polar and nonpolar interactions were color‐coded and indicated at the right bottom, which was built using Discovery Studio 2019. 2D, two dimensional

MD simulation

In computer‐aided drug design, MD simulation studies are used as a validation tool that predicts the binding stability of the complexes throughout the simulation run. Herein, the top three docked complexes with ARHGEF1 were examined via a 50 ns MD simulation run. Trajectories were extracted to analyze the dynamic properties of each receptor–ligand complex. Additionally, MD simulations of the top two ligands from COVID‐19 related library (Mitoxantrone and Camostat) complexed with ARHGEF1 were also performed for validation purposes. Each ligand shares the common binding pocket and exhibited similar kinds of interactions throughout the MD run. Additionally, five properties were also extracted from 50 ns MD trajectory for each hit compound, which includes (a) energy, (b) RMSD, (c) RMSF, (d) radius of gyration, and (e) hydrogen bonds.

The average potential energy score was found higher in ARHGEF1‐Lead3 as −7.91041e5 kcal/mole in comparison to control protein (ARHGEF1) which has an average potential energy of −8.76562e5 kcal/mole. The ARHGEF1‐Lead1 and ARHGEF1‐Lead2 had average potential energy of −1.16112e6 and −1.16114e6 which represents the stability of the complexes in terms of energy (Figure 5A). The potential energy of ARHGEF1‐mitoxantrone and ARHGEF1‐camostat complex were also energetically stable (Figure S7A). To measure the structural and conformational stability of the protein, the RMSD of the backbone atoms (Figure 5B and Figure S7B) for each complex was computed with respect to the initial docking structure. In all the receptor–ligand complexes and control protein, backbone RMSD showed deviation less than 0.4 nm. ARHGEF1‐Lead1, ARHGEF1‐Lead3, and ARHGEF1‐mitoxantrone showed similar kinds of trends like control protein, whereas ARHGEF1‐camostat and ARHGEF1‐Lead2 complex was the least stable complex in comparison to other complexes and control proteins. Furthermore, the flexibility of the amino acid residues was calculated using the RMSF profile. From Figure 5C and Figure S7C, RMSF of the backbone atoms was examined where residues of all complexes were found in the acceptable range of RMSF values, except in the case of F3394‐0943 against ARHGEF1 in which there is a minimal amount of fluctuation from Glu586 to Leu599, but not in the binding region of the protein. From Figure 5D and Figure S7D, it was observed that the average R g of all the complexes was below 2 nm which indicates the compactness of the system. Moreover, the number of hydrogen bonds that had a crucial role in the stability in the protein–ligand complex was analyzed from simulation trajectories. All the complexes remained stable (Figure 5E and Figure S7E), throughout the 50 ns run signified rigidity in the structure. The post‐MD binding interactions of the top three leads were analyzed, where all the hit compounds remain in the binding pocket of the protein (Figure 6).

Figure 5

Plots to investigate the energy conformations, stability, and fluctuations of ARHGEF1 apo and ARHGEF1 bound state: (A) represents energies, (B) represents RMSD, (C) represents RMS fluctuations, (D) radius of gyration, and (E) hydrogen bond. Black color shows the apo form of protein, red color, green color, and blue color represents Lead 1 (F3222‐4380), Lead 2 (F3394‐0943), and Lead 3 (F6548‐1087) bound ARHGEF1. RMS, root mean square; RMSD, root‐mean‐square deviation

Figure 6

Binding modes of the screened ligand into the active site of protein after 50 ns of MD run. DH, Dbl homology; MD, molecular dynamics

MM‐PBSA analysis

The MM‐PBSA binding free energy calculation of the top five complexes (including two‐hit compounds from COVID‐19 related library) were performed to revalidate the MD simulation and molecular docking exercises and are tabulated in Table 3. All the top five complexes exhibited negative binding energy in the MM‐PBSA threshold, respectively. However, among three compounds, Lead 1 (F3222‐4380) exhibits a better net binding energy score which may be comparable with the flexible CDOCKER energy score (−50.0296 kcal/mole) (Table 2). The difference in net binding energy among all the complexes is due to the difference in polar solvation energy score which plays a minimum role in protein–ligand interactions.

Table 3

MM‐PBSA energetic profile for the top three complexes

Compound Name	∆EVDW (kcal/mol)	∆EELE (kcal/mol)	∆GSOL‐PB (kcal/mol)	∆GSOL‐NP (kcal/mol)	∆Gbind‐PBSA (kcal/mol)
Lead 1 (F3222‐4380)	−50.35	−7.20	20.76	−4.13	−40.92
Lead 2 (F3394‐0943)	−48.89	−22.75	42.69	−4.47	−33.42
Lead 3 (F6548‐1087)	−45.51	−31.60	54.63	−4.02	−26.50

Abbreviations: MM, module membership; PBSA, Poisson–Boltzmann surface area.

DISCUSSION

COVID‐19 pandemic has radically affected the human population. Now, more than a year into the COVID‐19 has passed but still, no specific therapeutic option is available for its treatment. This is one of the main reasons that new cases and deaths are still reported throughout the world. Therefore, identification of new therapeutic targets and development of therapeutic agents in COVID‐19 are vigorously being pursued. On the basis of the life cycle of SARS‐CoV‐2 and the interaction of other molecular factors in COVID‐19, we can propose several targets for drug development as well as propose therapeutic agents (existing drug molecules and new ones) for the treatment of COVID‐19. This can be achieved through computational approaches that significantly speed up the identification of therapeutic targets and the development of therapeutic molecules in drug discovery. Many new therapeutic targets and therapeutic molecules have been proposed in COVID‐19. , , (p. 2), ,  But only a few have shown inhibition potency in bioassays, warranting new drug targets and drugs for the treatment of COVID‐19. Some of the repurposed drug molecules such as remdesivir have been registered for COVID‐19 therapy although based on the overall response of this drug globally till now states that it has several side effects. In addition to this, several vaccines have been approved for clinical use in the prevention of COVID‐19.

To serve the same purpose, we applied the WGCNA to identify biologically relevant therapeutic targets in PBMC COVID‐19 patient samples which were categorized into two groups, that is, Group 1: 5 mild + 5 controls and Group 2: 5 severe + 5 controls. WGCNA is a well‐known approach for gene co‐expression network analysis. , We specifically aimed to identify the common therapeutic target (i.e., MSI‐DEGs) in both groups. The DEGs overlapping between two groups was considered as MSI‐DEGs. We found 569 MSI‐DEGs which along with different age groups of patients were subjected to WGCN construction and hub module identification. In the WGCNA, we identified five unique gene modules. The green module (containing 56 MSI‐DEGs) which was significant and highly correlated (Tables S1–S3) was selected as our hub module. Among these MSI‐DEGs three upregulated (PPY, GNG4, and ARHGEF1) and one downregulated (CCL4L2) participated in the PPI network. While analyzing these MSI‐DEGs for the most significant pathway, GO‐BP, and GO‐MF terms, we found that 12 MSI‐DEGs participated in the top 10 significant pathways, 14 MSI‐DEGs participated in the top 10 GO‐BP, and 10 MSI‐DEGs participated in top 10 GO‐MF terms. The MSI‐DEG(s) overlapping between constructed highest confidence PPI network, top 10 significant pathways, and GO terms were considered as SARS‐CoV‐2 specific hub MSI‐DEG(s), respectively. We revealed ARHGEF1 as the overlapping hub MSI‐DEG between PPI, pathway, and GO terms (i.e., GO‐BP and GO‐MF). ARHGEF1 gene encodes for the RhoGEF1 protein. , RhoGEF1 is a member of a group of four RhoGEF proteins known to be activated by G protein‐coupled receptors (GPCRs) coupled to the G12 and G13 heterotrimeric G proteins. , We report that the most significant pathways associated with ARHGEF1 are G alpha (12/13) signaling events and signaling by nerve growth factor. The GO‐MF terms showed that ARHGEF1 is associated with Rho guanyl‐nucleotide exchange factor activity (GO: 0005089), while GO‐BP terms revealed its association with Rho protein signal transduction (GO:0007266), regulation of intracellular signal transduction (GO:1902531), and regulation of small GTPase mediated signal transduction (GO:0051056).

The involvement and overlapping of ARHGEF1 between PPI, significant pathways, and GO terms reflect its translational product as a suitable target for the therapeutics and therapeutic management of COVID‐19. Studies have reported that ARHGEF1, an intracellular signaling molecule, is involved in restraining signaling events from GPCRs to RhoA. It is crucially implicated in GPCR signaling events by possessing two regulatory domains – regulator of G‐protein signaling (RGS) domain that combines with Gα13‐containing heterotrimeric G‐proteins and mitigates GPCR signaling; and tandem DH/PH domains that participate in RhoA activation by encoding Rho guanine nucleotide exchange factor (GEF). Activation of RhoA, in turn, regulates several cellular processes and hence ARHGEF1 also regulates cell migration. (p. 13), (p. 13) Several reports have demonstrated its role in cell migration and motility. , Moreover, it also plays a critical role in pulmonary immunity. It was demonstrated by a study where ARHGEF1 existed as a tetramer in a B cell line, however, interrupting intramolecular association, lead to constitutive activation of RGS and RhoGEF. On the other hand, deficiency of ARHGEF1 increased adhesion of B lymphocytes, while its constitutive activation caused increased adhesion of B lymphocytes. In addition, it has also been found that ARHGEF1 is a key player in the production of Th2 cytokines by T cells, in the case of airway hyperresponsiveness and inflammation. Furthermore, ARHGEF1 deficient mice exhibited increased leukocytes and pulmonary components such as increased airspace and matrix metalloproteinase (MMP) activity with decreased lung elastase; typically found in patients with chronic obstructive pulmonary disease. , , , , It is also found to be involved in a novel signaling pathway, comprising of thromboxane A2 that is endowed in pulmonary inflammation. High or low levels of ARHGEF1 determines thromboxane receptor‐dependent production of MMP9 and macrophage association of fibronectin.

Additionally, there is no ligand complexed crystal structure available for ARHGEF1 at present. Likewise, there is no clear‐cut information available regarding the ligand‐binding sites present on ARHGEF1. Therefore, to target the ARHGEF1 for the therapeutics of COVID‐19, we analyzed its structure and predicted three binding pockets. Pocket 1 was defined as a suitable binding site based on rank score and composition of amino acids namely Thr427, Ala430, His431, Met434, Arg551, Leu552, Asp556, Met557, Thr560, Gln563, Arg564, and Lys567. The majority of these amino acid residues provide a plinth for the binding of potent inhibitors of many therapeutic targets in CoV diseases. (p. 2), ,

Next, we applied molecular docking and MD method against ARHGEF1 gene to discover potential hits for SARS‐CoV‐2 from the antiviral library, anticancer library, and natural product‐based library. We found that the three top ligands – Lead F3222‐4380, F3394‐0943, and F6548‐1087 possessed significant interactions with the active site of ARHGEF1, including hydrogen bonds, pi‐cation, and Van der Waals Interactions, suggesting their potential in modulating functional roles of ARHGEF1. To validate our findings and endorse the molecular docking study, we have selected the top 3 hit compounds for MD simulation. The RMSD profile from F3222‐4380, F3394‐0943, and F6548‐1087 did not exceed 0.4 nm as their averages were 0.2134, 0.2967, and 0.2336 nm, respectively, which represents their structural stability. The F3222‐4380 and F6548‐1087 complex were energetically more stable than the control protein. On the other hand, RMSF analysis exhibited that F3394‐0943 fluctuated more in the region of Glu586 and Leu599 in comparison to F3222‐4380 and F6548‐1087, overall, all the hit compounds have RMSF profile have less than 1 nm. These potential candidate hits could act as probable compounds that may target the important gene ARHGEF1, implicated in COVID‐19, and can be further selected for in vitro and in vivo studies in future work to speed up the drug discovery against SARS‐CoV‐2. Additionally, we performed molecular docking and dynamics method against COVID‐19 related library which serves as a positive control in our study. However, our novel hit compounds had better cDocker energy and MM‐PBSA profile in comparison to mitoxantrone and camostat (COVID‐19‐related compounds). Overall, our study has identified some important components involved in pulmonary immunity, and understanding the MMs behind the ARHGEF1 mediated regulation of immunity may reveal potential therapeutic targets to limit chronic inflammation. Besides, chemical compounds targeting ARHGEF1 could pave the way towards the development of medicine for the calamitous disease.

In conclusion, as we navigate through the findings, our study has identified some important components involved in pulmonary immunity. This study showed a significant association of the ARHGEF1 gene with COVID‐19, therefore highlighting it as a probable drug target for conventional treatments of this pandemic disease. The understanding of the MMs behind the ARHGEF1 mediated regulation of immunity may limit chronic inflammation in COVID‐19 patients. This study also found some actionable drug molecules with a rationale against the ARHGEF1 gene. These molecules may be effective alone or in combination with other drug molecules in COVID‐19. The in vitro and in vivo validation are warranted to determine their therapeutic efficacy against the devastating virus.

CONFLICT OF INTERESTS

The authors declare that there are no conflict of interests.

AUTHOR CONTRIBUTIONS

Prakash Jha, Prithvi Singh, Ravins Dohare, and Madhu Chopra: Conceptualization. Prakash Jha, Prithvi Singh, Arnab Nayek, and Kalaiarasan Ponnusamy: Methodology. Prakash Jha, Prithvi Singh, and Arnab Nayek: Software. Prakash Jha and Prithvi Singh: Formal analysis. Prakash Jha, Prithvi Singh, Arnab Nayek, and Kalaiarasan Ponnusamy: Data curation. Prakash Jha, Prithvi Singh, Armiya Sultan, and Shweta Arora: Writing – original draft. Prakash Jha, Prithvi Singh, Armiya Sultan, Shweta Arora, Arnab Nayek, Kalaiarasan Ponnusamy, Ravins Dohare, Mansoor Ali Syed, and Madhu Chopra: Writing – review and editing. Ravins Dohare and Madhu Chopra: Supervision. Ravins Dohare and Madhu Chopra: Project administration.

Supporting information.

Click here for additional data file.