Identification of Novel Gene Signatures using Next-Generation Sequencing Data from COVID-19 Infection Models: Focus on Neuro-COVID and Potential Therapeutics

SARS-CoV-2 is the causative agent for coronavirus disease-19 (COVID-19) and belongs to the family Coronaviridae that causes sickness varying from the common cold to more severe illnesses such as severe acute respiratory syndrome, sudden stroke, neurological complications (Neuro-COVID), multiple organ failure, and mortality in some patients. The gene expression profiles of COVID-19 infection models can be used to decipher potential therapeutics for COVID-19 and related pathologies, such as Neuro-COVID. Here, we used the raw RNA-seq reads (Single-End) in quadruplicates derived using Illumina Next Seq 500 from SARS-CoV-infected primary human bronchial epithelium (NHBE) and mock-treated NHBE cells obtained from the Gene Expression Omnibus (GEO) (GSE147507), and the quality control (QC) was evaluated using the CLC Genomics Workbench 20.0 (Qiagen, United States) before the RNA-seq analysis using BioJupies web tool and iPathwayGuide for gene ontologies (GO), pathways, upstream regulator genes, small molecules, and natural products. Additionally, single-cell transcriptomics data (GSE163005) of meta clusters of immune cells from the cerebrospinal fluid (CSF), such as T-cells/natural killer cells (NK) (TcMeta), dendritic cells (DCMeta), and monocytes/granulocyte (monoMeta) cell types for comparison, namely, Neuro-COVID versus idiopathic intracranial hypertension (IIH), were analyzed using iPathwayGuide. L1000 fireworks display (L1000FWD) and L1000 characteristic direction signature search engine (L1000 CDS2) web tools were used to uncover the small molecules that could potentially reverse the COVID-19 and Neuro-COVID-associated gene signatures. We uncovered small molecules such as camptothecin, importazole, and withaferin A, which can potentially reverse COVID-19 associated gene signatures. In addition, withaferin A, trichostatin A, narciclasine, camptothecin, and JQ1 have the potential to reverse Neuro-COVID gene signatures. Furthermore, the gene set enrichment analysis (GSEA) preranked method and Metascape web tool were used to decipher and annotate the gene signatures that were potentially reversed by these small molecules. In conclusion, our study unravels a rapid approach for applying next-generation knowledge discovery (NGKD) platforms to discover small molecules with therapeutic potential against COVID-19 and its related disease pathologies.


INTRODUCTION
Coronaviruses (CoVs) belong to the order Nidovirales, family Coronaviridae, and subfamily Coronavirinae, which can further be divided into four genera: alpha, beta, gamma, and delta CoVs. SARS CoV2 is the causative agent of coronavirus disease-19 , belongs to the genus beta-CoV, and can cause sickness varying from the common cold to more severe illnesses such as severe acute respiratory syndrome, gastrointestinal complications, sudden stroke, multiple organ failure, and mortality in some cases (Cui et al., 2019). SARS-CoV-2 infected more than 186 million people, resulting in the death of about 4 million people globally (Johns Hopkins COVID-19 Data Center on 10th July 2021) (Dong et al., 2020). SARS-CoV-2 has a positive-sense RNA genome encapsulated by a nucleocapsid. SARS-CoV-2 infects host cells through surface receptors, angiotensin-converting enzyme 2 (hACE2), and transmembrane protease serine-type 2 (TMPRSS2) (Hoffmann et al., 2020). An increase in the expression of ACE2, a tissueprotective mediator during lung damage, was found to be associated with interferon signaling in airway epithelial cells, and SARS-CoV-2 could exploit interferon-mediated stimulation of ACE2 to augment infection (Ziegler et al., 2020). The differential expression of genes that are necessary for SARS-CoV-2 interaction and subsequent host response determine susceptibility to COVID-19, disease progression, and recovery (Kasela et al., 2021).
RNA sequencing is a recently developed NGS methodology for whole transcriptome or single-cell transcriptomic approaches (Liu and Di, 2020). Single-cell RNA sequencing of COVID-19 infected bronchial epithelial cells and bronchioalveolar immune cells revealed important cellular and molecular processes implicated in COVID-19 infection at the single-cell level and provided information about the mechanisms of disease severity (Liu T. et al., 2020;Liao et al., 2020;Zhou et al., 2020). Notably, IL-17-associated signaling was significantly increased but not Th2-related inflammation following COVID-19 infection (Kasela et al., 2021). A recent study showed that SARS-CoV-2 infection caused a twofold higher induction of interferon stimulation compared to SARS-CoV in Calu-3 human epithelial cells and subsequent induction of cytokines such as IL6 or IL-10 (Wyler et al., 2021). The interferon-induced genes IFIT2 and OAS2 were widely stimulated compared to interferon lambda (IFNL) and interferon-beta (IFNB). Besides, scRNA-seq data suggested that interferon regulatory factor (IRF) activity occurs before the induction of nuclear factor-κB (NF-κB) in SARS-CoV-2infected cells (Wyler et al., 2021).
COVID-19 patients, especially those with greater disease severity, can develop neurological complications such as neuroinflammation, headache, and cerebrovascular disease called Neuro-COVID (Heming et al., 2021). Developing novel drug candidates and identifying suitable existing therapeutics for drug repurposing for COVID-19 and Neuro-COVID are critical for controlling this ongoing pandemic and reducing the enormous economic burden on health care systems and socioeconomic devastation of individuals, families, small to large businesses, and countries. Understanding COVID-19associated gene signatures is essential for developing robust therapeutics for treating infected patients effectively and reducing infection rates and mortality. To address this important issue, the gene expression profiles of COVID-19 infection models can be used to identify potential therapeutic targets that could be targeted by known drugs. Here, we used RNA-seq datasets from the COVID-19 infection model of human bronchial epithelial cells (NHBE) and the scRNA-seq datasets of immune cells isolated from the cerebrospinal fluid (CSF) of Neuro-COVID patients, obtained from public repositories and analyzed using next-generation knowledge discovery (NGKD) platforms to understand disease-specific gene signatures and uncover drugs from synthetic and natural sources that can reverse these gene signatures for potential therapeutics.

Ethical Statement
This study was exempted from Institutional Review Board (IRB) approval since it did not involve any animal models or human subjects and was conducted using RNA-seq datasets retrieved from the Gene Expression Omnibus (GEO) (Barrett et al., 2013).

Data Source
In the present study, the raw RNA-seq reads (Single-End) (FASTQ format) in quadruplicates derived using Illumina Next Seq 500 from SARS-CoV-infected and mock-treated NHBE cells were obtained from the GEO (Accession No: GSE147507) (Blanco-Melo, et al., 2020). Additionally, the single-cell transcriptomics data (Accession No: GSE163005) of immune cells isolated from the CSF of Neuro-COVID patients (Heming et al., 2021) were used for additional analysis using high-throughput knowledge discovery platforms. Heming et al. (2021) provided the entire dataset for the open-source interactive platform cerebroApp at http://covid.mheming.de/ (Hillje et al., 2020).

BioJupies Analysis of the RNA-seq Data
BioJupies is a freely available web-based application (http:// biojupies.cloud) that has 14 RNA-seq analysis library plug-ins and provides the user with the automatic generation, storage, and deployment of Jupyter Notebooks containing RNA-seq data analyses (Torre et al., 2018). In BioJupies, the RNA-seq datasets were user-submitted, compressed in an HDF5 data package, and uploaded to Google Cloud. Raw counts were normalized to log10-Counts per million (log CPM) and the differentially expressed genes (DEGs) were derived between the control group and the experimental group using the limma R package (Ritchie et al., 2015). The Jupyter Notebooks created for each RNA-seq raw data analysis were permanently available through a URL and stored in the cloud. The notebooks consist of executable code of the whole pipeline, description of the methods, enrichment analysis, interactive data visualizations, differential expression, and so on (Torre et al., 2018).
Principal component analysis (PCA) was performed using the PCA function from the sklearn Python module by transforming the log CPM using the Z-score method. An interactive heatmap was generated using a clustergram (Fernandez et al., 2017). In the volcano plots, the log2 fold changes of the DEGs are shown on the x-axis and p-values were corrected using the , and presented on the y-axis (Benjamini and Hochberg, 1995;Benjamini and Yekutieli, 2001). In contrast, for the MA plot, average gene expression is shown on the x-axis; p-values were corrected using the Benjamini-Hochberg method (Benjamini and Hochberg, 1995;Benjamini and Yekutieli, 2001), transformed (-log10), and presented on the y-axis.
In Silico Analysis of the RNA-seq Expression Data Using iPathwayGuide The impact analysis method (IAM) Khatri et al., 2007;Tarca et al., 2009) was used to determine the significantly impacted gene signatures and pathways from the DEGs (log2FC cut-off 0.6, adjusted false discovery rate (FDR) p-value ≤ 0.05) obtained from the COVID-19 using BioJupies and the DEGs with log2FC (cut-off 0.3) and adjusted p-value ≤ 0.001 based on Bonferroni method in meta clusters of T-cells/natural killer cells (NK) (TcMeta), dendritic cells (DCMeta), and monocytes/granulocyte (monoMeta) cell types of the comparison, namely, Neuro-COVID versus idiopathic intracranial hypertension (IIH) for the Neuro-COVID infection models in the iPathwayGuide (Advaita Bioinformatics, United States). Here, the p-value calculated based on Fisher's method was used to compute the pathway score method (Fisher, 1925). The p-value was further corrected based on multiple testing corrections for FDR and Bonferroni corrections (Bonferroni, 1935;Bonferroni, 1936). The gene interactions and pathways based on the DEGs were generated using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (Kanehisa and Goto, 2000;Kanehisa et al., 2002;Kanehisa et al., 2010;Kanehisa et al., 2012;Kanehisa et al., 2014). For each gene ontology (GO) term (Ashburner et al., 2000;Gene Ontology Consortium, 2001;Ashburner and Lewis, 2002;Gene Ontology Consortium, 2004), the number of DEGs annotated to the term was compared to that expected by chance. iPathwayGuide uses an overrepresentation approach to compute the statistical significance of observing at least a given number of DEGs (Draghici et al., 2003a;Draghici et al., 2003b;Draghici 2011). The hypergeometric distribution was used to compute the p-values in the iPathwayGuide analysis and corrected using FDR and Bonferroni for multiple comparisons (Draghici et al., 2003a;Draghici et al., 2003b;Draghici 2011). The prediction of upstream chemicals, drugs, and toxins (CDTs), either as present (or overly abundant) or absent (or insufficient), is based on two types of information: 1) the enrichment of DEGs from the experiment and 2) a network of interactions from the Advaita Knowledge Base (AKB v2012) (Draghici et al., 2003a;Draghici 2011). The analysis uses Fisher's standard method to combine p-values into one test statistic (Fisher, 1925).

L1000CDS2 and L1000FWD Queries
The L1000 characteristic direction signature search engine (L1000CDS2) analysis was performed by submitting the top 2000 DEGs to the L1000CDS2 signature search application programming interface (API) (Duan et al., 2016). Similarly, the L1000FWD analysis was performed by submitting the top 2000 DEGs to the L1000 Fire Works Display (L1000FWD) signature search API (Wang et al., 2018). Similarly, the DEGs obtained from TcMeta, DCMeta, and monoMeta cell types were compared; namely, Neuro-COVID versus IIH were subjected to both L1000CDS2 and L1000FWD analyses to identify drugs that reverse the gene signatures differentially regulated by COVID-19. An FDR (q-value) of 0.05 was considered statistically significant. desktop application GSEA 4.1.0 (Broad Institute, United States) under default settings as described previously (Subramanian et al., 2005).

Metascape Analysis of Gene Signatures Reversed by Small Molecules
The Metascape web tool (http://metascape.org) offers an easy and effective way to explore and understand gene lists derived from experimental data. The gene signatures reversed by small molecules identified in our study in COVID-19 and Neuro-COVID models were first automatically converted into Human Entrez Gene ID in Metascape. Then, all statistically enriched terms, accumulative hypergeometric p-values, and enrichment factors were calculated and used for filtering to obtain enrichment ontology clusters based on GO/KEGG terms, canonical pathways, and hallmark gene sets (Zhou et al., 2019).

RESULTS
Raw RNA-seq reads (Single-End) (FASTQ format) derived using Illumina Next Seq 500 from SARS-CoV-infected NHBE and mock-treated NHBE cells were obtained from the GEO and the QC was evaluated using the CLC Genomics Workbench 20.0, before the RNA-seq analysis using BioJupies web tool. iPathwayGuide analysis was performed to decipher the disease-specific signatures, pathways, and small molecules, either synthetic or derived from natural sources, to reverse disease-specific gene signatures. In addition, single-cell transcriptomic data of immune cells isolated from the CSF of Neuro-COVID-19 patients were further analyzed using iPathwayGuide, L1000CDS2, and L1000FWD analyses.
Hierarchically clustered heatmaps were generated using the Clustergrammer web tool to visualize and analyze highdimensional RNA-seq data of SARS-CoV-infected NHBE cells and mock-treated NHBE cells (Supplementary Figure S1A). PCA was used to uncover global patterns in RNA-seq datasets analyzed and helped to understand the difference between COVID-19-infected and mock-treated NHBE cells (Supplementary Figure S1B). The volcano plot was generated using transformed gene fold changes using log2 and is shown on the x-axis (Supplementary Figure S1C). The MA plot was based on the average gene expression, which was calculated using the mean of the normalized gene expression values and shown on the x-axis (Supplementary Figure S1D).

iPathwayGuide Analysis of DEGs From COVID-19 and Neuro-COVID Infection Models
In this experiment, 1,072 DEGs were identified from a total of 10,663 DEGs obtained from BioJupies analysis of the RNA-seq reads of the SARS-CoV-infected NHBE cells and mock-treated NHBE cells based on a p-value cut-off of 0.05 and a log2 fold change cut-off of 0.6. In contrast, DEGs with a logFC cut-off of 0.3 and adjusted p-value based on the Bonferroni method from clusters in TcMeta, DCMeta, and monoMeta of the comparison, namely, Neuro-COVID versus IIH, were also subjected to iPathwayGuide analysis separately, followed by comparative analyses. Subsequently, the DEGs were analyzed in the context of pathways obtained from the KEGG database (Kanehisa and Goto, 2000;Kanehisa et al., 2002), GO from the Gene Ontology Consortium database, a network of regulatory relationships from BioGRID: Biological General Repository for Interaction Datasets v4.0.189 (Szklarczyk et al., 2017), chemicals/drugs/toxicants from the Comparative Toxicogenomics Database (Davis et al., 2021), and diseases from the KEGG database. In summary, 22 pathways were found to be significantly impacted in SARS-CoV-2-infected NHBE cells compared to mock-treated NHBE cells. In addition, 503 GO terms, 18 miRNAs, 190 gene upstream regulators, 213 chemical upstream regulators, and 14 diseases were found to be significantly enriched before correction for multiple comparisons.
COVID-19 infection of NHBE cells triggers key immunerelated pathways, such as cytokine-cytokine receptor interactions and viral protein interactions with cytokine receptors ( Table 1). The top five upstream regulators, IL-17, TNF-alpha, STAT2, IRF9, and TLR4, were predicted to be activated ( Table 2). The top identified biological processes, molecular functions, and cellular components for each pruning type are provided in Tables 3-5.
The bar chart ( Figure 1A) shows the top small molecules identified by the L1000CDS2 query using the DEGs identified from SARS-CoV-2-NHBE. The left panel shows small molecules such as geldanamycin, radicicol, AZD8330, trametinib, NVP-AYU922, GSK2126458, and JW-7-24-1, which mimic the observed gene expression signature; the right panel displays small molecules such as camptothecin ( Figure 1B), importazole ( Figure 1C), and withaferin A ( Figure 1D). The upstream regulator drugs and natural products that reverse the molecular signatures based on iPathwayGuide analysis are shown as a dendrogram ( Figure 1E). The top five upstream drugs, natural products, and chemicals predicted as absent (or insufficient) based on iPathwayGuide analysis were coumestrol, methylprednisolone, JinFuKang (JFK), selenium, and gold sodium thiomalate ( Figure 1F). However, withaferin A was found to reverse the COVID-19-induced molecular signatures in both L1000CDS2 and L1000FWD analyses, along with other small-molecule drugs ( Table 6) in the SARS-CoV-2-infected NHBE cells.  Table 2 indicates the number of differentially expressed (DE) targets supporting the hypothesis that each upstream regulator (u) is activated (DTA(u)), the total number of DE genes downstream of u (DT(u)), the combined raw p-value, and the corrected p-value for multiple comparisons. We identified 14 genes that were commonly expressed between Neuro-COVID and IIH (TcMeta, DCMeta, and monoMeta), as depicted in the Venn diagram ( Figure 2A). The upregulated genes ( Figure 2B), downregulated genes ( Figure 2C), and the common genes between the meta clusters of immune cells in Neuro-COVID were presented as rank diagrams based on log2FC values. The genes GABARAP, GNAI2, COTL1, ATP5F1D, CD81, GNAS,   gas transport, oxygen transport, hydrogen peroxide, and drug transport ( Figure 3B), the top five molecular functions enriched were haptoglobin binding, oxygen binding, oxygen carrier activity, heme binding, and tetrapyrrole binding ( Figure 3D), and the top five cellular components enriched were endocytic vesicle, haptoglobin-hemoglobin complex, Frontiers in Pharmacology | www.frontiersin.org August 2021 | Volume 12 | Article 688227 8 hemoglobin complex, cytosolic small ribosomes, and cytosolic ribosomes ( Figure 3F). The top five differentially expressed pathways identified based on iPathwayGuide analysis of immune cell meta clusters from Neuro-COVID patients were malaria, African trypanosomiasis, cocaine addiction, Parkinson's disease, and leukocyte transendothelial migration. The differentially regulated pathways in the meta clusters of immune cells from patients with Neuro-COVID are provided in Supplementary Table S1. The upstream genes activated in TcMeta, DCMeta, and monoMeta clusters are listed in Supplementary Table S2.

L1000FWD and L1000CDS2 Analyses
The L1000FWD analysis of DEGs of meta clusters of Tc, DC, and Mono of Neuro-COVID compared to IIH revealed that withaferin A was the top molecule capable of reversing the COVID-19 induced gene signatures (Tables 7-9). Furthermore, the rank diagram ( Figure 4A) showed that JQ1 was the top drug based on the in silico prediction of insufficient signaling of drugs, natural products, and chemicals in the meta clusters of Tc, DC, and Mono of Neuro-COVID compared to IIH using iPathwayGuide ( Figure 4B). The L100CDS2 analysis of DEGs of meta  clusters of Tc, DC, and Mono of Neuro-COVID compared to IIH revealed that narciclasine ( Figure 4C) and trichostatin A ( Figure 4D) were some of the top molecules potentially reversing the Neuro-COVID gene signatures (Tables  10-12).

DISCUSSION
COVID-19 caused by SARS-CoV-2 infection remains an ongoing pandemic (Huang C. et al., 2020;Liu J. et al., 2020;Novel Coronavirus Pneumonia Emergency Response Epidemiology Team, 2020) and patients with severe COVID-19 may also develop neurological complications called Neuro-COVID (Heming et al., 2021). RNA sequencing is a very recently developed NGS methodology for the whole transcriptome or single-cell transcriptomics approaches (Liu and Di, 2020) and is broadly used to explore biological, cellular, and molecular processes implicated in COVID-19 infection (Liu T. et al., 2020;Liao et al., 2020;Zhou et al., 2020). Hence, either developing novel drug candidates or identifying suitable existing therapeutics for drug repurposing for COVID-19 and Neuro-COVID is essential to decrease the infection rate and control the COVID-19 pandemic and reduce the enormous economic burden on healthcare systems. Because the gene expression profiles of COVID-19 infection models can be used to decipher potential therapeutic targets that could be targeted by known drugs (Daamen et al., 2021), we used RNA-seq datasets from the COVID-19 infection models of NHBE cells, and the scRNA-seq datasets of immune cells isolated from the CSF of Neuro-COVID patients and analyzed using NGKD platforms to understand the disease-specific gene signatures and pathways and further uncover small molecules from both synthetic and natural sources that potentially reverse these diseases.
Here, we found that COVID-19 infection of NHBE cells activated upstream genes such as IL-17, TNF-alpha, STAT2, IRF9, and TLR-4. Biological processes such as humoral immune response, acute-phase response, and molecular functions such as cytokine activity, receptor regulator activity, signaling receptor activity, receptor-ligand activity, and chemokine activity were enriched in the COVID-infected cells. Importantly, the cytokine and cytokine receptor interaction and viral protein interaction with the cytokine and cytokine receptors were activated in COVID-infected NHBE cells. Cytokines are important for both innate and adaptive inflammatory host responses, cell differentiation, cell death, growth, repair and development, and cellular homeostasis (Pushparaj, 2019;Bahlas et al., 2020;Harakeh et al., 2020;Jafri et al., 2020;Pushparaj, 2020). Studies have shown that several circulating cytokines and chemokines such as TNFα, CXCL-10, IL-6, and IL-8 are differentially expressed during SARS-CoV-2 infection, and this cytokine/chemokine storm likely contributes to the poor prognosis of COVID-19 (Liu J. et al., 2020;Vaninov, 2020). RNA sequencing analysis of cell and animal models of SARS-CoV-2 infection, blood, lung, and airway biopsies from COVID-19 patients showed inflammatory responses characterized by low levels of type I and III IFNs, increased interleukin-6 (IL-6), and a variety of chemokines (Blanco-Melo et al., 2020;Daamen et al., 2021). The spike protein (S protein) of SARS-CoV-2 is essential for the attachment between the coronavirus and hACE2 surface receptor through its receptor-binding domain (RBD) (Lan et al., 2020) and is proteolytically activated by human proteases, thus helping the coronavirus to enter the host cells (Shang et al., 2020). A recent study showed that hACE2 was stimulated by IFN in human airway epithelial cells (Ziegler et al., 2020) and thus helps in the entry of SARS-CoV-2 into host cells.
SARS-CoV-2 and other coronaviruses have developed different mechanisms to avoid detection and subsequent destruction by copying and repurposing cytokine and cytokine receptor genes in the host (Heimfarth et al., 2020;Choudhary et al., 2021). COVID-19 induced cytokines and cytokine receptors, chemokines, and other specific cytokine receptors and binding proteins may subvert and alter the host cytokine networks (Choudhary et al., 2021). Here, the COVID-19-induced cytokines, cytokine receptors, receptor-binding proteins, and chemokines may stimulate or prevent cytokine signaling and may significantly alter various facets of host immunity. In addition, Daamen et al. (2021) found that COVID-19 pathogenesis was driven by highly inflammatory myeloidlineage cells with distinct transcriptional signatures and the absence of cytotoxic cells in the lungs, leading to reduced viral clearance. Heming et al. (2021) stated that lumbar puncture to obtain immune cells from COVID-19 patients without neurological manifestations as controls was not ethically permitted for scientific purposes. Since IIH is a benign disorder associated with high pressure in the brain, the immune cells derived from the CSF of patients with IIH were used as controls to compare Neuro-COVID. The cluster of differentiation molecule 81 (CD81) is one of the commonly regulated genes in the meta clusters of immune cells from the CSF of patients with Neuro-COVID and belongs to the tetraspanin superfamily, which has been shown to regulate viral entry, viral replication, infectivity, and virion exit of different types of viruses (Benayas et al., 2020). Therefore, it is essential to investigate the importance of CD81 in patients with COVID-19 and Neuro-COVID. One of the upstream genes activated in TcMeta cluster, Cell cycle division 37 (CDC37), a heat shock protein 90 (HSP90) cochaperone that could play an important role in the pathogenesis of Neuro-COVID. COVID-19 progression to a systemic disease could be associated with HSP-related molecular mimicry autoimmune phenomena (Cappello et al. 2020;Kasperkiewicz, 2021). It was postulated that Hsp90 inhibition could also be a potential treatment option for cytokine storm-mediated acute respiratory distress syndrome in COVID-19 patients (Kasperkiewicz, 2021). Recently, Wyler et al. (2021) identified HSP90 as a target for COVID therapy based on transcriptomic profiling of SARS-CoV-2 infected human cell lines.
Interestingly, the top five differentially expressed pathways identified based on iPathwayGuide analysis of immune cell meta clusters from Neuro-COVID patients were malaria, African trypanosomiasis, cocaine addiction, Parkinson's disease, and leukocyte transendothelial migration. Studies have shown a potential link between the presentation of malaria and COVID-19. The opposite relationship between COVID 19 and malaria has been suggested to be linked with the wide use of antimalarial drugs, including hydroxychloroquine (HCQ) and chloroquine (CQ), in countries that are endemic to malaria (Hussein et al., 2020).
There are many types of COVID-19 vaccines currently available for prophylaxis, and many are under development (Mandolesi et al., 2021). Several therapeutics are available based on WHO guidelines to treat the complications of COVID-19 and related complications (Lamontagne et al., 2021); however, these therapeutics are not specifically designed for the treatment of COVID-19 and its related complications such as Neuro-COVID, and their efficacies substantially differ across the globe and are not very effective in ameliorating disease severity (Surnar et al., 2020). In this study, we utilized NGKD platforms such as iPathwayGuide, L1000FWD, and L1000CDS2 tools to identify promising druggable molecules based on their in silico potential to reverse gene signatures induced by COVID-19 and Neuro-COVID. We found that camptothecin, importazole, and withaferin A had insufficient signaling or gene signatures (or absent) in COVID-19 infected NHBE cells. Based on L1000CDS2 analysis, trichostatin A, a histone deacetylase inhibitor, mildly inhibited the ACE receptors (Takahashi et al., 2021), and narciclasine and camptothecin are some of the top small molecules that reverse the gene signatures in Neuro-COVID vs. IIH immune datasets. In addition, a comparative analysis of the Neuro-COVID vs. IIH immune cell meta cluster datasets showed that JQ1 had insufficient signaling (or absence).
The GSEA preranked analysis calculates if a priori defined sets of genes display statistically significant enrichment at either end of the ranking (Subramanian et al., 2005). The gene signature potentially reversed by withaferin A in SARS-CoV-2 NHBE vs. Mock-NHBE based on preranked GSEA involved in various biological, molecular, and cellular processes, including viral genome replication (GO: 0019079), modulation of the process of other organisms involved in symbiotic interactions (GO:0051817), and positive regulation of translational initiation (GO:0045948). The gene signature potentially reversed by importazole in SARS-CoV-2 NHBE vs. Mock-NHBE based on preranked GSEA involved in various biological, molecular, and cellular processes, including regulation of single-stranded viral RNA replication via a double-stranded DNA intermediate (GO: 0045091). The gene signature potentially reversed by withaferin A in Neuro-COVID vs. IIH-TcMeta based on preranked GSEA involved in various biological, molecular, and cellular processes, including regulation of type 1 interferon production (GO:0032479) and interferon signaling (R-HSA-913531). The gene signature potentially reversed by withaferin A in Neuro-COVID vs. IIH-TcMeta based on preranked GSEA involved in various biological, molecular, and cellular processes, including regulation of type 1 interferon production (GO:0032479) and interferon signaling (R-HSA-913531). The gene signature potentially reversed by narciclasine in Neuro-COVID vs. IIH-TcMeta based on preranked GSEA involved in negative regulation of viral entry into host cells (GO: 0046597), PDGFR beta signaling pathway (PID-M186), EGF/ EGFR signaling pathway (WP437), VEGFA/VEGFR2 signaling pathway (WP3888), and positive regulation of cell migration (GO: 0030335). The gene signatures upregulated by withaferin A (GNLY CST7 PPIB TSPO BCL2 S100A10 GSTP1), (S100A10, TSPO, PPIB, HLA-DQB1, BCL2, GSTP1, EDF1, and FLOT1), and S100A9, S100A8, TSPO, and HOMER3) were positively enriched in Neuro-COVID. Granulysin (GNLY) is a member of the saposinlike protein (SAPLIP) family, is located in the cytotoxic granules of T-cells and NK cells, is released on antigen stimuli, and has antimicrobial activity. The S100 genes include 13 members and have antibacterial and antifungal properties (Crinier et al., 2018).
Our L1000FWD analyses showed that withaferin A was the top natural product that reverses the signature of Neuro-COVID in all the meta clusters of immune cells from the CSF of Neuro-COVID patients. Withaferin A is a component of Withania somnifera (ashwagandha or Indian ginseng) (Srivastava et al., 2020). W. somnifera has been used in traditional medicine as an antioxidant, antianxiety, anti-inflammatory, antibacterial, aphrodisiac, and herbal tonic for general health (Sood et al., 2018). The active ingredients include withanolides, saponins, alkaloids, and steroidal lactones. In vitro studies have shown that ashwagandha has neuroprotective, cardioprotective, immunomodulating, and anticancer properties .
Adjunctive treatment with ashwagandha improved symptoms and stress in patients with schizophrenia, offering beneficial effects on cognitive function in patients with bipolar disorder and improves balance in patients with progressive degenerative cerebral ataxias (Sood et al., 2018;Singh et al., 2021). It was recently shown that withanolides present in ashwagandha possess anti-COVID-19 properties, and these compounds exhibit good absorption and transport kinetics with no related mutagenic or adverse effects (Srivastava et al., 2020). Withaferin A was predicted to bind and stably interact with the binding site of TMPRSS2, similar to its known inhibitor, camostat mesylate (Kumar et al., 2020). Camostat was found to reduce SARS-CoV-2 infection in TMPRSS2 expressing Vero cells (Hoffmann et al., 2020). David et al., 2021 in their MedRixv preprint showed that a common variant of TMPRSS2 protects against COVID-19. In silico screening of several phytochemicals identified that Withanone, one of the constituents of ashwagandha, showed a potential inhibition of ACE2 (Balkrishna et al., 2021). Additionally, Ghosh et al. (2021) used molecular dynamic simulations and pharmacophore modeling approaches to predict the highly potent small-molecule derivative of withaferin A that potentially inhibits SARS-CoV-2 protease (Mpro), a favorable future therapeutic against  Recent studies have demonstrated the antiviral properties of narciclasine, an alkaloid found in various Amaryllidaceae species, and camptothecin, a topoisomerase inhibitor first isolated from the stem of Camptotheca acuminata (used in Chinese traditional medicine) against SARS-CoV-2 (Huang C.-T. et al., 2020;Mamkulathil Devasia et al., 2021). However, importazole, an inhibitor of importin-β transport receptors, and other small molecules identified to reverse COVID-19-induced gene signatures need to be further explored because developing effective therapeutics is essential to control the COVID-19 pandemic (Surnar et al., 2020).
In conclusion, the present study unravels a rapid approach to using high-throughput RNA sequencing technologies coupled with NGKD platforms to decipher specific drugs and small molecules derived either synthetically or from natural sources for the amelioration of COVID-19 related disease pathologies such as Neuro-COVID. Further studies are warranted to validate the small molecules identified in our study using in vitro and in vivo model systems of COVID-19 and Neuro-COVID to determine their mechanism(s) of action followed by suitable clinical trials to confirm the efficacy and safety for possible therapeutic intervention for COVID-19related disease pathologies.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found here: Gene Expression Omnibus GSE147507 and GSE163005.

AUTHOR CONTRIBUTIONS
PP, MN, and AA designed the experiments. PP and MN conducted the experiments. PP, AA, and MN analyzed the data. PP wrote the manuscript. PP and MN finally revised the manuscript. All authors contributed to the editing of the manuscript and the scientific discussions.

FUNDING
The research work was funded by Institutional Fund Project under grant no. (IFPHI-110-117-2020). Therefore, the authors gratefully acknowledge technical and financial support from the Ministry of Education and King Abdulaziz University, DSR, Jeddah. Saudi Arabia.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphar.2021.688227/ full#supplementary-material Supplementary Table S1 | The differentially regulated pathways in the meta clusters of immune cells from patients with Neuro-COVID are provided in Supplementary Table S1.
Supplementary Table S2 | The upstream genes activated in TcMeta, DCMeta, and monoMeta clusters are listed in Supplementary Table S2.
Supplementary Datasheet S1 | The genes enriched in GSEA preranked analysis of SARS CoV2-NHBE vs. Mock-NHBE against the gene signatures of camptothecin, importazole, and withaferin A are provided in Supplementary Datasheet S1.
Supplementary Datasheet S2 | The genes enriched in GSEA preranked analysis of Neuro-COVID vs. IIH comparison against the gene signatures of withaferin A, camptothecin, trichostatin A, narciclasine, and JQ1 are provided in Supplementary Datasheet S2.
Supplementary Figure S1 | (A) Hierarchically clustered interactive heatmaps were generated using the Clustergrammer web tool for visualizing and analyzing highdimensional RNASeq data (NHBE-SARS CoV2 vs NHBE-Mock). (B) Principal Component Analysis (PCA) was applied to identify global patterns in highdimensional RNASeq datasets (C) Volcano plot was generated using transformed gene fold changes using log2 and displayed on the x-axis (D) MA plot was based on average gene expression which was calculated using the mean of the normalized gene expression values and displayed on the x-axis.
Supplementary Figure S6 | GSEA Preranked analysis was performed to decipher the potential gene signatures differentially regulated (Combined, Downregulated, and Upregulated) by narciclasine using the RNK file generated from DEGs of Neuro-COVID vs IIH (TcMeta, DCMeta, and monoMeta) comparisons. The enrichment of gene signature (Signature ID: Frontiers in Pharmacology | www.frontiersin.org August 2021 | Volume 12 | Article 688227