CSF3 Is a Potential Drug Target for the Treatment of COVID-19

Coronavirus Disease 2019 (COVID-19) is an acute respiratory infectious disease that appeared at the end of 2019. As of July 2020, the cumulative number of infections and deaths have exceeded 15 million and 630,000, respectively. And new cases are increasing. There are still many difficulties surrounding research on the mechanism and development of therapeutic vaccines. It is urgent to explore the pathogenic mechanism of viruses to help prevent and treat COVID-19. In our study, we downloaded two datasets related to COVID-19 (GSE150819 and GSE147507). By analyzing the high-throughput expression matrix of uninfected human bronchial organoids and infected human bronchial organoids in the GSE150819, 456 differentially expressed genes (DEGs) were identified, which were mainly enriched in the cytokine–cytokine receptor interaction pathway and so on. We also constructed the protein–protein interaction (PPI) network of DEGs to identify the hub genes. Then we analyzed GSE147507, which contained lung adenocarcinoma cell lines (A549 and Calu3) and the primary bronchial epithelial cell line (NHBE), obtaining 799, 460, and 46 DEGs, respectively. The results showed that in human bronchial organoids, A549, Calu3, and NHBE samples infected with SARS-CoV-2, only one upregulated gene CSF3 was identified. Interestingly, CSF3 is one of the hub genes we previously screened in GSE150819, suggesting that CSF3 may be a potential drug target. Further, we screened potential drugs targeting CSF3 by MOE; the top 50 drugs were screened by flexible docking and rigid docking, with 37 intersections. Two antiviral drugs (Elbasvir and Ritonavir) were included; Elbasvir and Ritonavir formed van der Waals (VDW) interactions with surrounding residues to bind with CSF3, and Elbasvir and Ritonavir significantly inhibited CSF3 protein expression.


INTRODUCTION
Coronavirus Disease 2019, first detected in Wuhan, China, in December 2019, has been raging around the world from its detection through to current times (Kannan et al., 2020;Thompson, 2020). On January 9, 2020, the World Health Organization (WHO) announced the discovery of the new coronavirus. This coronavirus was originally called 2019-nCoV and then officially named SARS-CoV-2, which had never been found among humans. On February 11, 2020, the respiratory disease originating from SARS-CoV-2 infection was called COVID-19 Mahase, 2020). Although the virus was first discovered in Wuhan, there is still controversy about its specific origin (Andersen et al., 2020;Forster et al., 2020). WHO announced that COVID-19 had become a global health problem because it has spread in more than 140 countries around the world (Zhai et al., 2020).
Fever, sore throat, fatigue, cough, or dyspnea are the typical symptoms of COVID-19 Ozma et al., 2020). As of 10 August 2020, confirmed COVID-19 cases have reached 19,718,030, including 728,013 deaths in the world 1 . The government shared the genome sequence of SARS-CoV-2 with the public in a timely manner. This was effective in helping scientists deal with the emergency (Harcourt et al., 2020;Wu et al., 2020). Although most countries in the world have effectively controlled COVID-19, it is still far from being completely over. This is because research on the treatment and prevention of COVID-19 is still insufficient, and the vaccine under research has not yet been fully applied in clinic (Dhama et al., 2020).
At present, the treatment of COVID-19 is still limited, and it has been reported that some drugs may have a therapeutic effect on COVID-19. As an antimalarial agent, chloroquine and its derivative, hydroxychloroquine, have shown to be active against COVID-19 in vitro (Wang M. et al., 2020). However, there is not enough clinical data to support the therapeutic effect of chloroquine or hydroxychloroquine on COVID-19 Huang et al., 2020).
In addition, remdesivi (Jean et al., 2020), lopinavir, and emetine (Choy et al., 2020) have been reported to inhibit SARS-CoV-2 in vitro. But these drugs are not currently recommended for clinical use, and are only recommended for use in clinical trials (Kotwani and Gandra, 2020). Moreover, research on the screening of anti-SARS-CoV-2 drugs and their mechanism is mainly focused on angiotensin-converting enzyme 2 (ACE2) (Williamson and Kerimi, 2020). Whether there are other drug targets for SARS-CoV-2 remains an unanswered question for researchers.
In our study, we analyzed the high-throughput sequencing data of uninfected human bronchial organoids and SARS-CoV2 infected human bronchial organoids in the GEO database to explore the effects of SARS-CoV2 (Blanco-Melo et al., 2020). Then we verified the high-throughput data of A549, Calu3, and NHBE cell lines, finding that only CSF3 (Tamada et al., 2006) was up-regulated in the four sets of samples. It has been reported that the granulocyte colony-stimulating factor (CSF3) was the most upregulated gene after SARS-CoV2-infected genes (Nunnari et al., 2020), suggesting that CSF3 is a significant target after infection and it may also be a potential target for drug therapy. Therefore, we have screened drugs that may potentially target CSF3 for the treatment of COVID-19 from the FDAapproved drug library.

Data Download
The Gene Expression Omnibus (GEO) 2 database was used to downloaded datasets about COVID-19. With "COVID-19 OR SARS-CoV-2, Homo sapiens" as the screening criteria, we selected GSE150819 and GSE147507 for analysis. In GSE150819, human bronchial organoids (hBO) from commercially available cryopreserved primary human bronchial epithelial cells (hBEpC) were constructed for severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) research. We chose three uninfected samples and three SARS-CoV-2 infected samples to analyze. In GSE147507, we selected primary human lung epithelium cells (NHBE) and lung adenocarcinoma cells (A549, Calu-3), which were mock treated or infected with SARS-CoV-2.

Identification of DEGs
In our study, | log fold change (FC)| > 1.5 and p < 0.05 were used as the cut-off criteria to find DEGs in uninfected and SARS-CoV-2-infected organoids or cell lines. Then, these DEGs were visualized by volcano plots and heatmaps, which were performed by R package in RStudio. VENNY 2.1.0 3 online software was used to identify overlapping DEGs in the four sets of samples.

GO and KEGG Enrichment Analysis
Gene Ontology (GO) analysis classified the genes in the difference tables according to their functions, which included biological process (BP), cellular component (CC), and molecular function (MF). The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis showed the key pathways of clustered genes. GO and KEGG enrichment analysis of DEGs in GSE150819 were performed by the Database for Annotation, Visualization, and Integrated Discovery (Huang et al., 2009) (DAVID) 4 v6.8. And a p < 0.05 was considered to be statistically significant.

Protein-Protein Interaction (PPI) Network Construction, Module Analysis, and Identification of Hub Genes
In order to facilitate understanding of the interaction between the proteins corresponding to these DEGs, a PPI network was constructed by Search Tool for the Retrieval of Interacting Gene (STRING) 5 database (Szklarczyk et al., 2015). In the network, core modules and top 10 hub genes were identified by ClusterONE and CytoHubba in Cytoscape software, suggesting the key genes and modules which might play a significant role in COVID-19.

Gene Set Enrichment Analysis
Gene set enrichment analysis (GSEA) is a software which can analyze the positive and negative regulatory pathway of a hub gene or all DEGs. To better understand the role of CSF3 in COVID-19, we chose KEGG function of CSF3 in GSE10819 to analyze, obtaining the key positive and negative regulatory pathways by NOM p-value ranking (p < 0.05).

Molecular Docking
The X-ray structure of CSF3 was downloaded from RCSB Protein Data Bank (PDB ID: 2D9Q) 6 . The FDA-approved drug library was obtained from DrugBank 7 . The 2D structures of the FDAapproved drug molecules were drawn in ChemBioDraw 2014 and converted to 3D in MOE through energy minimization. MOE dock was used for homology modeling and molecular docking simulations of the binding affinity between FDA-approved drugs and the target protein. Prior to docking, the force field of AMBER10:EHT and the implicit solvation model of Reaction Field (R-field) were selected. The protonation state of the protein and the orientation of the hydrogens were optimized by LigX, at a PH of 7 and temperature of 300 K. The docking workflow followed the "induced fit" protocol, in which the side chains of the receptor pocket were allowed to move according to ligand conformations, with a constraint on their positions. The weight used for tethering side chain atoms to their original positions was 10. All docked poses were ranked by London dG scoring first, then a force field refinement was carried out on the top 10 poses followed by a rescoring of GBVI/WSA dG. The conformation with the lowest free energy of binding was selected as the best (probable) binding mode. Molecular graphics were generated by PyMOL.

Western Blotting
Crude cellular proteins were separated by SDS-PAGE and transferred to nitrocellulose membranes. Membranes were blocked in 5% non-fat milk and probed with primary antibodies (1:1,000 dilution) against CSF3 (Abcam, ab181053) overnight at 4 • C. Membranes were incubated with horseradish peroxidase (HRP)-conjugated secondary antibodies for 1.5 h and immunolabeling detected using a Bio-Rad imaging system.

Identification of DEGs in Datasets
In GSE150819, we analyzed three samples of uninfected human bronchial organoids and three samples of SARS-CoV-2-infected human bronchial organoids, obtaining 456 DEGs, of which 205 were upregulated DEGs and 251 downregulated DEGs (Table 1).
To visualize the expression of the DEGs in GSE150819, a heatmap and volcano plot were constructed, as shown in Figure 1.

Enrichment Analysis of DEGs
To identify the significant pathways associated with the DEGs between uninfected human bronchial organoids and infected human bronchial organoids in GSE150819, we performed the GO and KEGG pathway enrichment analysis by DAVID. The result showed that in GO enrichment analysis ( Table 2), genes were enriched in biological processes (BP) of inflammatory response, immune response, and cell chemotaxis. As for cellular component (CC), these genes showed enrichment in extracellular space, extracellular region, and integral component of plasma membrane. The molecular function (MF) results showed cytokine activity, chemokine activity, and structural molecule activity. KEGG pathway enrichment analysis was shown in Table 2. DEGs were mainly enriched in cytokine receptor interaction, TNF signaling pathway, and African trypanosomiasis.

Construction of PPI Network and Identification of Hub Genes
STRING online software constructed the PPI network of 456 DEGs, composed of 282 nodes and 930 edges after excluding the isolated nodes, as shown in Figure 2A. ClusterONE identified the top three modules in the PPI network, suggesting that these submodules might be highly significant in the process of COVID-19. Module 1 ( Figure 2B) and module 2 (Figure 2C), consisting of 22 and 12 DEGs, respectively, were both enriched in KEGG pathway, including the Cytokine-cytokine receptor interaction and TNF signaling pathway. Module 3 ( Figure 2D) contained 14 nodes and 45 edges.

Identification of Overlapping Genes in Different Datasets and GSEA Analysis
Similarly, we also analyzed the three cell lines samples, which were NHBE, A549, and Calu3, and their SARS-CoV-2-infected cells in GSE147507. After screening, the results suggested that there were 45 DEGs in the NHBE group, 158 DEGs in the A549 group, and 460 DEGs in the Calu3 group. VENNY 2.1 online software was used to find the overlapping DEGs in the four sets of samples, and finally only one overlapping DEG was found, as shown in Figures 3A,B. Only CSF3 showed a common upregulated gene in infected human bronchial organoids and cell lines. For GSEA, we analyzed CSF3 enrichment in GSE150819 to explore its potential molecular functions, obtaining the positive and negative regulatory pathways of CSF3. The results were shown in Figure 3C. The results showed that the top four positive regulatory pathways of CSF3 were proteasome, Parkinson's disease, oxidative phosphorylation, and Graft vs. host disease. While the top four negative regulatory pathways of CSF3 were mismatch repair, DNA replication, homologous recombination, and inositol phosphate metabolism.

Using CSF3 as the Target to Find Potential Targeted Drugs
Giuseppe Nunnari discovered the potential value of CSF3 in COVID-19 research (Nunnari et al., 2020). We obtained the X-ray structure of CSF3 (Tamada et al., 2006) from the RCSB PDB and aimed to screen out the 2,470 FDA-approved drugs to potentially target CSF3. We investigated the binding mode of FDA-approved drugs with CSF3 by rigid and flexible docking simulation. The scores of the top 50 in flexible docking are shown in Table 3. Among them, ubiquinol (DB11340), an active antioxidant (Ernster and Forsmark-Andrée, 1993;Sies, 1997), got the highest score of -8.0513706 kcal/mol. In addition, two antiviral drugs, Elbasvir (DB11574) (Myers et al., 2015) with scores of -6.5620656 kcal/mol and Ritonavir (DB00503) (Cameron et al., 1998) with scores of -6.259393 kcal/mol, have attracted our attention.
The binding mode of Elbasvir with CSF3 were illustrated in Figure 4A. Elbasvir has formed a suitable steric complementarity with the binding site of CSF3. The oxygen atom of carbonyl in Elbasvir, regarded as hydrogen bond acceptors, form the hydrogen bonds with Lys 16 and Lys 23 in CSF3. The binding mode of Ritonavir with CSF3 were illustrated in Figure 4B. The nitrogen atom in Ritonavir, regarded as a hydrogen bond donor, forms hydrogen bonds with Gln 119 in CSF3. Moreover, Elbasvir and Ritonavir formed van der Waals (VDW) interactions with FIGURE 3 | Identification and analysis of overlapping DEGs. VENNY software identified the upregulated DEGs (A) and downregulated DEGs (B) of the four sets of samples in GSE150819 and GSE147507. (C) The overlapping gene CSF3 was analyzed by GSEA. There were four positive regulatory pathways of CSF3, which were proteasome, Parkinson's disease, oxidative phosphorylation, and Graft vs. host disease. The following were four negative regulatory pathways of CSF3: mismatch repair, DNA replication, homologous recombination, and inositol phosphate metabolism.
surrounding residues, which mainly contribute to the binding energy between Elbasvir and Ritonavir with CSF3.
We further verified the effect of Elbasvir and Ritonavir on CSF3; WB results show that both Elbasvir and Ritonavir significantly inhibited the CSF3 protein expression in A549 cells (Figures 4C,D).

DISCUSSION
In the current global pandemic, COVID-19 has infected more than 19 million people and caused more than 720,000 deaths. It has a very high infection rate , and the case fatality rate has been reported to be as high as 5% (Li L.Q. et al., 2020). Clinical data in the United States showed that among patients hospitalized in the intensive care unit, case fatality was up to 40% (Wiersinga et al., 2020). The virus has also bene shown to spread among people of all ages 8 . Studies have shown that the risk of death is highly correlated with the age gradient, and younger age may be associated with lower mortality (Verity et al., 2020).
This is an epidemic that has affected all mankind and requires our joint efforts. In the winter when the virus was raging (Smit et al., 2020), China quickly quarantined Wuhan, promptly and effectively suppressing the spread of COVID-19. The Chinese people came together with the government to combat the disease.
As of March 2020 9 , the daily number of local cases in mainland China has dropped to two digits, and the infected are mainly from people returning from abroad. But the virus still exists, and we still have not completely defeated it. The existing obstacles are the limitations of treatment methods, the lack of effective drugs (Triggle et al., 2020), the difficulty of vaccine development (Uddin et al., 2020), and the inability to completely block the spread of the virus (Rothan and Byrareddy, 2020). All of these issues require continuous efforts by researchers, and research into treatment and prevention should be given equal attention. On the one hand, reliable treatments and therapeutic drugs need to be found. On the other hand, the spread of the virus needs to be blocked and vaccines need to be researched as quickly as possible.
In our study, we hope to find potential drug therapeutic targets from the high-throughput data of virus-infected organoids and cells and explore whether there are clinical drugs in use that can target the infectious protein. We found 456 DEGs in the analysis of uninfected and SARS-CoV-2-infected human bronchial organoids. KEGG enrichment analysis showed that these genes were mainly enriched in Cytokine-cytokine receptor interaction and TNF signaling pathway. Through the verification of high-throughput data of A549, Calu3, and NHBE cell lines, we have found an important target protein CSF3, which, along with its receptor CSF3R, participate in the regulation of granulopoiesis, neutrophil function, and hematopoietic stem cell mobilization (Zhang et al., 2018). It has been reported that CSF3 mutations have carcinogenic effects (Hollmén et al., 2016;Jing et al., 2016). Virus-related research has found that human respiratory syncytial virus (hRSV) infection disrupted the polarity of the pediatric respiratory epithelial secretome and was associated with immune modulating proteins (CXCL6, CXCL16, and CSF3) never before linked with this virus before (Touzelet et al., 2020). The latest reports showed that CSF3 is the most modulated gene in NHBE cells infected with SARS-CoV-2 (Nunnari et al., 2020), indicating the potential value of CSF3.
Regarding CSF3 as a drug target, we use molecular docking simulation technology by MOE to screen out potential drugs from the FDA-approved drug library, which were downloaded from DrugBank. These drugs were scored separately by rigid docking and flexible docking. Among the top 50 drugs, a total of 37 drugs were screened by both methods. The screened drugs included anti-HCV drugs, such as Ombitasvir (DB09296) (Ascione et al., 2018), Daclatasvir (DB09102) (Lee, 2013), and Elbasvir (DB11574) (Myers et al., 2015), and the anti-HIV drugs, Cobicistat (DB09065) (Angione et al., 2018) and Ritonavir (DB00503) (Cameron et al., 1998). The drug with the highest score is Ubiquinol (DB11340), which is an active antioxidant drug (Ernster and Forsmark-Andrée, 1993). No antiviral drugs have been proven to be effective against COVID-19. Therefore, we give priority to antiviral drugs such as Elbasvir to explore whether it has a potential effect to COVID-19 in our study. Similarly, Junmei Wang screened out Elbasvir, carfilzomib, eravacycline, valrubicin, and lopinavir as potential inhibitors of SARS-CoV-2's main protease (Wang, 2020). In Meenakshisundaram Balasubramaniam's study, Elbasvir was predicted to bind multiple SARS-CoV-2 proteins, blocking virus replication, and had a good antiviral effect (Balasubramaniam and Reis, 2020). As an FDA-approved drug for anti-HIV and HCV, Ritonavir (Feld et al., 2014) was tried early in the treatment of COVID-19 (Lim et al., 2020;Wang Y. et al., 2020), although the results were not satisfactory (Griffin, 2020). At present, these predicted antiviral drugs are still not recommended for clinical use, and are only supported for clinical trials. In a word, further mechanism exploration and clinical trials are needed for the application of these antiviral drugs.

CONCLUSION
In conclusion, our study identified the DEGs and underlying pathways in uninfected and SARS-CoV-2-infected human bronchial organoids and cell lines, found the only upregulated gene was CSF3, and explored its potential value. We regarded CSF3 as a potential drug target, and used molecular docking and virtual drug screening technology to find drugs that may interact with it to find potential treatments for COVID-19. Elbasvir, Ritonavir, and other antiviral drugs were screened through both rigid docking and flexible docking, and they were potential anti-COVID-19 drugs that have been studied. Targeting CSF3 may be a potential therapeutic mechanism of these drugs, but a more in-depth exploration is needed.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
CF, NW, WZ, and Y-LL: conception and design. NW and QL: administrative support. JM, CF, and HT: collection and assembly of data. CF, DR, and JM: data analysis and interpretation. NW, JM, and CF: manuscript writing. All authors finally approved the manuscript.