key: cord-323737-6ajqy0ch
authors: Jiang, Yuanyuan; Liu, Lanxin; Manning, Morenci; Bonahoom, Madison; Lotvola, Aaron; Yang, Zhe; Yang, Zeng-Quan
title: Structural analysis, virtual screening and molecular simulation to identify potential inhibitors targeting 2'-O-ribose methyltransferase of SARS-CoV-2 coronavirus
date: 2020-10-04
journal: Journal of biomolecular structure & dynamics
DOI: 10.1080/07391102.2020.1828172
sha: 
doc_id: 323737
cord_uid: 6ajqy0ch

SARS-CoV-2, an emerging coronavirus, has spread rapidly around the world, resulting in over ten million cases and more than half a million deaths as of July 1, 2020. Effective treatments and vaccines for SARS-CoV-2 infection do not currently exist. Previous studies demonstrated that nonstructural protein 16 (nsp16) of coronavirus is an S-adenosyl methionine (SAM)-dependent 2’-O-methyltransferase (2’-O-MTase) that has an important role in viral replication and prevents recognition by the host innate immune system. In the present study, we employed structural analysis, virtual screening, and molecular simulation approaches to identify clinically investigated and approved drugs which can act as promising inhibitors against nsp16 2′-O-MTase of SARS-CoV-2. Comparative analysis of primary amino acid sequences and crystal structures of seven human CoVs defined the key residues for nsp16 2-O’-MTase functions. Virtual screening and docking analysis ranked the potential inhibitors of nsp16 from more than 4,500 clinically investigated and approved drugs. Furthermore, molecular dynamics simulations were carried out on eight top candidates, including Hesperidin, Rimegepant, Gs-9667, and Sonedenoson, to calculate various structural parameters and understand the dynamic behavior of the drug-protein complexes. Our studies provided the foundation to further test and repurpose these candidate drugs experimentally and/or clinically for COVID-19 treatment. Communicated by Ramaswamy H. Sarma

Three coronaviruses (CoVs): severe acute respiratory syndrome coronavirus (SARS-CoV-1), Middle East respiratory syndrome coronavirus (MERS-CoV), and the recently identified SARS-CoV-2 in December 2019, have caused deadly pneumonia in humans since the beginning of the 21st century (Coronaviridae Study Group of the International Committee on Taxonomy of, 2020; Wu, Zhao, et al., 2020; Zhu et al., 2020) . The SARS-CoV-2 causes coronavirus disease-19 (COVID-19) with influenza-like symptoms ranging from mild discomfort to severe lung injury and multi-organ failure, eventually leading to death. As of July 1, 2020, more than ten million COVID-19 cases were reported worldwide, and more than half million patients have died (https://www.who. int/emergencies/diseases/novel-coronavirus-2019). Effective treatments and vaccines for SARS-CoV-2 infection do not currently exist. Thus, there is an urgent and unmet need for the discovery and development of antiviral therapeutics for the treatment of COVID-19.

CoVs are positive-sense RNA viruses that replicate in the cytoplasm of infected cells. Replication and transcription of the CoV RNA genome are achieved by a complex RNA replication/transcription machinery, consisting of at least 16 viral nonstructural proteins (nsps) (Coronaviridae Study Group of the International Committee on Taxonomy of, 2020; Perlman & Netland, 2009) . Previous studies demonstrated that nsp16 proteins of SARS-CoV-1 and MERS-CoV have methyltransferase (MTase) activities that catalyze methylation of the first transcribed nucleotide at the ribose 2 0 -O position (2 0 -O-Me) (Aouadi et al., 2017; Chen et al., 2011; Decroly et al., 2011) . The 2 0 -O-Me of viral cap RNA protects itself from degradation by 5 0 -3 0 exoribonucleases, ensures efficient translation, and helps to prevent recognition by the host innate immune system . The importance of 2 0 -O-MTase activity for CoV infection and pathogenesis was previously documented by in vitro and in vivo studies (Gonzales-van Horn & Sarnow, 2017; Sevajol et al., 2014; Subissi et al., 2014) . For SARS-CoV-1, the absence of nsp16 2 0 -O-MTase activity results in significant attenuation characterized by decreased viral replication, reduced weight loss, and limited breathing dysfunction in mice . In addition, 2 0 -O-Me acts as a recognition marker that helps the host cell recognize its own RNA species (Byszewska et al., 2014; Z€ ust et al., 2011) . Inhibition of nsp16 2 0 -O-MTase activities should restrain viral replication and enable recognition by the host innate immune system. This makes the nsp16-MTase a promising target for the identification of new anti-SARS-CoV-2 drugs.

Recent advances in structural bioinformatics and virtual screening approaches have revolutionized the identification and/or repurposing of marketed drugs or bioactive compounds for effective treatment of various human diseases, including infectious diseases (Chang et al., 2016; Gioia et al., 2017; Kitchen et al., 2004; Lasko et al., 2017; Maia et al., 2020; Pinzi & Rastelli, 2019; Slater & Kontoyianni, 2019) . Moreover, in silico approaches, including molecular dynamics (MD), have also been widely used to determine the conformational space of the investigated targets, ligands, and ligandtarget complexes, and thus better understand the dynamic behavior of ligand-target complexes Pinzi & Rastelli, 2019; Rajendran et al., 2018; Slater & Kontoyianni, 2019; Sled z & Caflisch, 2018) . Very recently, a virtual screening approach was used to identify potential drugs to inhibit SARS-CoV-2 proteins, including surface spike glycoprotein, main protease, and nsp16 Panda et al., 2020; Tazikeh-Lemeski et al., 2020; Vijayan et al., 2020) . For example, Bhardwaj et al. reported the identification of bioactive molecules from tea plant as SARS-CoV-2 main protease inhibitors . Molecular docking and virtual screening approaches also have been attempted to identify compounds targeting 2 0 -O-MTase of SARS-CoV-2 (Encinar & Menendez, 2020; Hijikata et al., 2020; Tazikeh-Lemeski et al., 2020; Vijayan et al., 2020) .

In the present study, we employed structural analysis, virtual screening, and molecular simulation approaches to identify potential inhibitors targeting 2 0 -O-MTase of SARS-CoV-2. We first performed comparative analysis of primary amino acid sequences and crystal structures of seven human CoVs and defined the key residues for nsp16 2 0 -O-MTase functions. We executed virtual screening and docking analysis to rank the potential inhibitors of nsp16 from more than 4,500 clinically investigated and approved drugs. MD simulations were carried out on eight candidate compounds to calculate various structural parameters and understand the dynamic behavior of the drug-protein complexes. Our studies provided the foundation to further test and repurpose these candidate drugs experimentally and/or clinically for COVID-19 treatment.

To identify inhibitors targeting nsp16, we first performed comparative analysis of primary amino acid sequences and crystal structures of seven human CoVs. Supplementary Table 1 lists the detailed genome and protein information that were employed in this study. In primary amino acid sequences, nsp16 of SARS-CoV-2 was found to be 93.3% identical to SARS-CoV-1, but only 56.6 À 65.9% identical to five other human CoVs (MERS-CoV, HCoV-OC43, -HKU1, -NL63, and -229E). The nsp16 proteins belong to a class of S-adenosyl methionine (SAM) -dependent 2 0 -O-MTases present in all life forms, and they all contain the conserved catalytic KDKE motif (K46, D130, K170, and E203 in SARS-CoV-2) ( Figure 1A & B) (Bouvet et al., 2010; Chen et al., 2011; Decroly et al., 2011) .

We next analyzed crystal structures of nsp16s obtained from SARS-CoV-2, SARS-CoV-1, and MERS-CoV that are available in RCSB Protein Data Bank (PDB) as of April 30, 2020 (Supplementary Table S1 ) (Berman et al., 2000; Chen et al., 2011; Decroly et al., 2011) . The three-dimensional (3 D) structures of nsp16s were solved in complex with nsp10 that functions as a stimulatory factor to regulate 2 0 -O-MTase function (Supplementary Figure S1 ) (Chen et al., 2011; Decroly et al., 2011) . Structural alignment was performed using TMalign, an algorithm for sequence-independent protein structure comparisons (Zhang & Skolnick, 2005) . We revealed the highest coverage, percent identity, and structural conservation between SARS-CoV-2 and SARS-CoV-1 nsp16 bound to natural ligand SAM [PDB: 6W4H chain A (Resolution: 1.80 Å); and 3R24 chain A (Resolution: 2.00 Å), respectively], with TMscore 0.958 and Root Mean Square Deviation (RMSD) ¼ 0.89. The SARS-CoV-2 nsp16 structure is also highly similar to the MERS-CoV complex bound to SAM and a cap analogue: 7methyl-GpppA [PDB:5YNM: chain A (Resolution: 1.68 Å)], with TM-score 0.956 and RMSD ¼ 0.80 (Supplementary Figure S2 ). The superposition of 6W4H and 5YNM revealed that the conserved KDKE motif of SARS-CoV-2 nsp16 is located at the bottom of the RNA binding pocket ( Figure 1C & D). Structural analysis also revealed that residues N43, Y47, G71, G81, D99, L100, N101, D114, and M131 coordinate SAM binding in both SARS-CoV-1 and À2 through hydrogen bonds and water-mediated interactions (Aouadi et al., 2017; Viswanathan et al., 2020) .

In humans, cap methyltransferase 1 (CMTR1) is a specific 2 0 -O-MTase that catalyzes methylation of the first transcribed nucleotide (N1) in cellular mRNA (Cap 1) (Belanger et al., 2010; Smietanski et al., 2014) . We also performed comparative analysis of the crystal structures of CoV-2 nsp16 and human 2 0 -O-MTase CMTR1. SARS-CoV-2 nsp16 and human CMTR1 MTase domain share little primary amino acid sequence identity (12.1%) based on the Pairwise Sequence Alignment with the Needleman-Wunsch algorithm at EMBL-EBI services (https://www.ebi.ac.uk/Tools/psa/). A search of the PDB using DALI service (http://ekhidna2.biocenter.helsinki.fi/dali/) revealed structural similarity between CoV-2 nsp16 (PDB: 6W4H chain A) and human CMTR1 (PDB: 4N48 and 4N49), which have Dali Z-Scores 13.7 and 13.8, respectively, and RMSD 3.4 (Supplementary Figure S3) . However, even though both nsp16 and CMTR1 contain a KDKE catalytic tetrad motif and a conserved folding made of a Rossmann b-sheet surrounded by helices on each side, substantial differences are found in the RNA binding groove as well as the SAM binding pocket ( Figures 1C and 2) . In CMTR1, Cap 0 binds in a deep pocket, with a triphosphate bridge and N1 phosphate backbone stabilized by a stack of arginine residues along the RNA binding groove (Smietanski et al., 2014) . In contrast, the RNA binding groove on SARS-CoV-2 nsp16 appears to be less positively charged and shallower ( Figure 1C ). Additionally, the SAM binding pocket on SARS-CoV-2 nsp16 is more negatively charged than that of human CMTR1. Also, the SAM binding pocket on SARS-CoV-2 nsp16 is more open and flexible compared to that of CMTR1, which is deeper and longer ( Figures 1C and 2) . The difference between the sequences and structures of the human 2 0 -O-MTase CMTR1 and CoV nsp16 could be exploited to develop inhibitors that could specifically block viral MTases.

We hypothesize that the small-molecule antagonists of nsp16 of SARS-CoV-2 will limit viral replication and/or unmask viral RNA to intracellular innate immunity. To identify new drugs that have the potential to inhibit nsp16 functions of SARS-CoV-2, we first examined the druggability of all binding sites of 2 0 -O-MTase by the DoGSiteScorer tool (Volkamer et al., 2012) . With values between zero and one, the SAM pocket has the highest drug score (0.85) and RNA binding pocket has a score of 0.55 in SARS-CoV-2 nsp16 (Supplementary Figure S1 ). Similar scores were observed in MERS-CoV nsp16 (PDB: 5YNM) that was solved in complex with both SAM and the cap RNA analogue (Supplementary Figure S1 ). As the aim of this work was to identify clinically investigated and approved drugs for quickly repurposing to treat COVID-19 infections, we made use of an MTiOpenScreen database Drugs-lib that contains 7,173 stereoisomers corresponding to 4,574 drugs, generated from ChEMBL, DrugBank, DrugCentral, and SuperDrug2 (Lagarde et al., 2018 (Lagarde et al., , 2019 . Among them, the nsp16 natural ligand SAM would serve as a positive control for virtual screening. Among the five available nsp16 structures of SARS-CoV-2 (April 30, 2020), we selected the PDB 6W4H structure which presents the best resolution (1.8 Å) to screen the Drugs-lib using the MTiOpenScreen service (Lagarde et al., 2019) . The computations of MTiOpenScreen were carried out with AutoDock Vina, a gradient optimization algorithm, which provided a list of the 1,500 highest-scoring compounds (Trott & Olson, 2010) . We repeated the MTiOpenScreen screening protocol three times using the 6W4H structure to ensure the reliability of the results. We found that four stereoisomers of SAM were ranked in the top 1,500 best compounds in all three runs, with energy equal to À8.10 kcal/mol ( Figure 3A & B, Supplementary Table  S2) . Among the 1,500 top hit compounds in each screening, 1,380 compounds were shared across all three MTiOpenScreen runs with highly similar AutoDock Vina scores. Based on the means of scores, we ranked these 1,380 compounds (stereoisomers) that correspond to 967 drugs (Supplementary Table  S2 ). The top ten drugs with the highest scores were MK3207, Rimegepant, Entrectinib, Osi-027, Bolazine, R428, Hesperidin, Losulazine, Rebastinib, and Cep-32496, with the means of energy around À10.97 $ À10.17 kcal/mol ( Figure 3A and Table 1) . Notably, there are 38 Hesperidin stereoisomers and four MK3207 stereoisomers among the 1,380 hits ( Figure 3A and Supplementary Table S2) .

To understand previously reported pharmacological action of hit compounds, we queried the biological target and pathway classes in the Probes and Drugs Portal, in which 713 standardized (693 non-isomeric) drugs were available for these 967 drugs (Skuta et al., 2017) . The top target classes of these drugs are G-protein-coupled receptors (241), catalytic receptors (139), and kinases (129). The top pathways are signal transduction (453), immune system (270), and transcription (247) ( Figure 3C ). For example, MK3207 and Rimegepant are calcitonin gene-related peptide (CGRP) receptor antagonists, which have been clinically investigated and/or approved for migraine treatment (Bell, 2014; Lipton et al., 2019; Salvatore et al., 2010) . Hesperidin possesses biological and pharmacological properties as an effective antioxidant, anti-inflammatory, anti-carcinogenic, and anti-hypertensive agent (Aggarwal et al., 2020; Tejada et al., 2018) .

Using the principles of structural similarity searching, analogues of the natural ligand (SAM) in nsp16 might be potent inhibitors against 2 0 -O-MTase function. Since SAM was among our 1,380 highest-scoring compounds, we next computed the similarity between each pair of compounds using the FragFp descriptors in DataWarrior (Sander et al., 2015) . Among them, we found that seven drugs are chemically similar to SAM, namely Gs-9667, Trabodenoson, Binodenoson, Sonedenoson, Regadenoson, Metrifudil, and Selodenoson ( Figure 3B and Table 2 ). All seven of these drugs are adenosine receptor agonists that were clinically investigated for treating various diseases, including cardiac arrhythmias, neuropathic pain, inflammatory diseases, and cancer (Gao & Jacobson, 2011; Jacobson et al., 2019) .

Rigid docking performed using Autodock Vina provided insight into the binding affinity and hydrogen bond formation of the candidate compounds on the nsp16 2 0 -O-MTase. Next, the top ten drugs and seven drugs structurally similar to SAM were chosen for the ligand flexible docking simulation by using the RosettaLigand program (Davis & Baker, 2009; Lemmon & Meiler, 2012) . RosettaLigand employs the Monte Carlo minimization protocol in which the ligand position and orientation are randomly perturbed by a small deviation (0.1 Å and 3 ). A scoring function in RosettaLigand includes an electrostatics model, an explicit orientation-dependent hydrogen bonding potential, an implicit solvation model, and van der Walls interactions (Davis & Baker, 2009; Lemmon & Meiler, 2012 ). Rosetta's energy scores (REU) of these 17 drugs were shown in Tables 1 and 2. With the RosettaLigand simulation, Hesperidin ranked the highest with the lowest REU among the top ten hits obtained from the Autodock Vina screening. In the seven SAM-like drugs, similar to AutoDock Vina results, Gs-9667 had the lowest REU score (Table 2) .

We next examined the interaction of these 17 drugs with key residues for nsp16 2 0 -O-MTase function. We found that Hesperidin, Osi-027, Gs-9667, and Sonedenoson have many hydrogen bonds or hydrophobic interactions with key functional residues involved in 2 0 -O-MTase function, which are represented in Figure 4 . Hesperidin forms hydrogen bonds with N43, K46, G71, L100, N101, C115, Y132, and K170 and hydrophobic interactions with Y132, and F149. Osi-027 forms hydrogen bonds with N43, K46, L100, N101, C115, Y132, and K170, a salt bridge with D99, and hydrophobic interactions with L100, D114, M131, and F149. Gs-9667 forms hydrogen bonds with L100, C115, and Y132 and salt bridges with D99 and D114. Sonedenoson forms hydrogen bonds with N43, K46, D75, N101, Y132, and K170, salt bridges with D75 and D99, and hydrophobic interactions with L100, M131, and F149. Notably, Hesperidin, Osi-027 and Sonedenoson form hydrogen bonds with two key residues in the KDKE motif: K46 and K170.

To gain insight into the dynamic behavior of the compounds at the active site of 2 0 -O-MTase, a set of 100 ns MD simulations were executed for eight selected hits (top hits and/or clinically approved drugs, Table 3 ) by employing the GROningen MAchine for Chemical Simulations (GROMACS) package (Van Der Spoel et al., 2005) . The obtained results were read as RMSD, which is essential to quantify the structural stability of protein-drug complexes within a regular time frame. RMSD analysis of the free form of 2 0 -O-MTase revealed that the protein molecule achieves stability at around 20 ns, and it mostly maintained its stability until 100 ns. Nsp16-SAM and nsp16-Hesperidin complexes attained stability at around 50 ns and maintained that stability until the end. Based on average RMSD, all eight drugs were very comparable to SAM, and the mean values of RMSD in all compound-protein complexes are lower than those in free protein (Table 3 ). The lowest mean value of RMSD is observed for the Gs-9667 simulation (Table 3) .

Root Mean Square Fluctuation (RMSF) is used to estimate the average atomic mobility of backbone atoms (N, Ca, and C) during MD simulation. The 2 0 -O-MTase nsp16 of SARS-CoV-2 has an extensively looped structure, and we observed the fluctuations throughout the entire nsp16 protein ( Figure 5 ). Even though the average fluctuations of all eight proteindrug complexes were very similar, the nsp16-Hersperidin complex was found to have the highest, while nsp16-Entrectinib complex had the lowest average RMSF value (Table 3) .

We also calculated additional structural parameters from the MD simulations, like Radius of Gyration (Rg), Solvent-Accessible Surface Area (SASA), Inter-molecular Hydrogen Bonding (H-bonding), and Principal Component Analysis (PCA) based on essential dynamics (ED) approach (Figures 5  and 6 ). Based on the Rg analysis that estimates the compactness of protein and protein-drug complexes, the 2 0 -O-MTase was undergoing some compression in its 3 D conformation ( Figure S4 ). Among the protein-drug complexes, the most compressed complex was that of Hesperidin (Table 3) . These results are consistent with the SASA analysis, which determined that the amount of overall surface area is reduced for solvent accessibility due to compression in the protein structure ( Figure 5) . Additionally, our MD simulation analysis showed an average of five hydrogen bonds present between nsp16 protein and two drugs (Hesperidin and Regadenoson), while others form an average of one to three hydrogen bonds ( Figure 5 and Table 3 ).

The PCA was employed to reveal relevant motions from the global trajectories of unbound proteins as well as the protein-drug complexes using the ED approach. ED projects the combined fluctuations of the most unsteady regions of protein molecules into two variables, namely Principal Component 1 (PC1) and Principal Component 2 (PC2). The projections of trajectories for free protein as well as complexes are illuminated in Figure 6 . Among them, the Hesperidin-protein complex was found to be occupying the most, while the Osi-027 complex had the least conformational space (Table 3) .

Next, we performed the MM-PBSA (Molecular Mechanics/ Poisson-Boltzmann Surface Area) analysis to estimate the protein-drug binding affinity in dynamic state. As shown in Table 4 , we found that three compounds, Hesperidin, Rimegepant, and Entrectinib, had the highest binding affinities with protein in dynamic state. Among three SAM-like drugs with MD simulation analysis, Sonedenoson and Gs-9667 had relatively higher binding affinities, compared to the SAM with 2 0 -O-MTase in dynamic state. To determine the energy contributions of individual residues to the total binding free energy, energy decomposition analysis was conducted on a per-residue basis and the result is presented in Supplementary Table S3 . As illuminated in Figure 7 , the calculations of residue-based free energy decomposition identified the hot interaction spots of 2 0 -O-MTase with four compounds. For example, four residues (D75, L100, Y132, and F149) of nsp16 were notable hotspots (< À10.0 kcal/ mol) that contribute favorably to the binding of Hesperidin.

In the Gs-9667 complex, two residues, D130 and Y132, of nsp16 were dramatically favorable toward binding (-10.6 and À20.2 kcal/mol). As seen from MM-PBSA results and docking studies, drugs including Hesperidin, Osi-027, Rimegepant, Sonedenoson, and Gs-9667 had higher binding affinities than SAM with the 2 0 -O-MTase of SARS-CoV-2.

Methylation of RNA, which occurs in all kingdoms of life, plays important roles in RNA metabolism, processing, stability, nuclear export, translation efficiency, and others (Kadumuri & Janga, 2018; Li & Mason, 2014; Roundtree et al., 2017; Shi et al., 2019) . 2 0 -O-Me is predominantly found in rRNA and tRNA of bacteria and eukaryotes, as well as in the 5 0 mRNA Cap of higher eukaryotes (Ayadi et al., 2019; Boccaletto et al., 2018) . In humans, five 2 0 -O-MTases [FtsJ RNA methyltransferase homolog 1 (FTSJ1), FTSJ2, FTSJ3, CMTR1, and CMTR2] have been reported to catalyze 2 0 -O-Me of various RNAs. These 2 0 -O-MTases regulate important cellular functions, and are associated with developmental disorders and cancer (B€ ugl et al., 2000; Hager et al., 2002; Lee & Bogenhagen, 2014; Ringeard et al., 2019) . Mutations of FTSJ1, a tRNA 2 0 -O-MTase, cause autosomal-recessive intellectual disability (Freude et al., 2004; Guy et al., 2015; . FTSJ2 (also known as MRM2) is a mitochondrial rRNA methyltransferase (Lee et al., 2013) . Human FTSJ3, which likely catalyzes 2 0 -O-Me of both rRNA and internal sites of mRNA, is a potential regulator of breast cancer progression (Manning et al., 2020) . Furthermore, FTSJ3 catalyzes 2 0 -O-Me of HIV RNAs and leads to the inhibition of innate immune sensing and response (Ringeard et al., 2019) . Both CMTR1 and CMTR2 catalyze mRNA cap methylation, CMTR1 for the first transcribed nucleotide and CMTR2 for the second transcribed nucleotide (Inesta-Vaquera & Cowling, 2017) . Notably, CMTR1 has previously been identified as ISG95, a protein implicated in the response to interferon treatment and viral infection (Geiss et al., 2003; Su et al., 2002) . A very recent study in the influenza A virus revealed that loss of CMTR1 in host cells inhibits viral replication and up-regulates anti-viral genes (Li, Clohisey, et al., 2020) . Thus, 2 0 -O-MTases have critical roles in viral replication and anti-viral immune response. Many RNA viruses, including CoVs, replicate in the cytoplasm using their own viral 2 0 -O-MTases, which catalyze the formation of cap structures on viral mRNA that mimic those present on host mRNAs (Hyde & Diamond, 2015; Netzband & Pager, 2020) . Previous studies on SARS-CoV-1 and MERS-CoV demonstrated that targeting nsp16 2 0 -O-MTases might interfere with viral replication by inhibiting the replication process and promoting intracellular recognition/immune response to viral RNA species Sevajol et al., 2014; Subissi et al., 2014) . Thus, nsp16 is a promising target for identification of drugs for viral infections, including SARS-CoV-2.

In this study, we performed virtual screening and ranked the predicted binding affinities of approximately 1,000 therapeutic drugs that have the potential to inhibit nsp16 2 0 -O-MTase function. We identified top ten candidates and seven SAM-like compounds. These drugs have been experimentally and clinically investigated for treatment of various diseases. For example, Hesperidin is one of the main components of a Chinese medicine, beniol (Bai et al., 2019) . Hesperidin had demonstrated antiviral properties in in vitro models against influenza A virus (H1N1) by stimulating cell-autonomous immunity through interferon gene expression (Ding et al., 2018) . Rimegepant recently received FDA approval for the acute treatment of migraine in adults (Scott, 2020) . Additionally, Entrectinib and Osi-027 are both anti-cancer drugs (Drilon, 2019; Mateo et al., 2016; Sartore-Bianchi et al., 2020) . SAM-like compound Gs-9667 was clinically investigated for treating Type 2 diabetes (Staehr et al., 2013) . To additionally compare the features of these candidate compounds in Tables 1 and 2, we also performed pharmacophore analysis using an online PharmaGist server (Schneidman-Duhovny et al., 2008) . The number of features and spatial feature set for each compound were summarized in Supplementary Table S4 . Multiple aligned structures of the top ten compounds revealed that the highest PharmaGist score is 29.7, which resulted from the alignment of six compounds (Supplementary Figure S5) . These six compounds shared four pharmacophoric features: three aromatic rings and one hydrogen bond acceptor. Furthermore, as expected, the seven SAM-like compounds have a relatively higher PharmaGist score. The highest score (50.6) resulted from the alignment of four compounds (Sonedenoson, Regadenoson, Metrifudil, Selodenoson) plus SAM. They shared ten features including two aromatic rings, two hydrogen bond donors, five hydrogen bond acceptors, and one positive charge (Supplementary Figure S6) . These four compounds and SAM likely form hydrogen bond with residue N101 of nsp16 protein. In the future, after determining how some of these compounds directly bind nsp16 of SARS-CoV-2 with various biochemical approaches, the pharmacophore model will help to identify highly potent compounds.

To gain insight into the dynamic behavior of the small molecules at the active site of 2 0 -O-MTase, we performed MD simulation analysis for eight top hits including MK3207, Hesperidin, Osi-027, Entrectinib, Rimegepant, Gs-9667, Sonedenoson, and Regadenoson. The results suggest that all eight compound-protein complexes were overall stable with acceptable deviation (Figure 5 ). Among the eight complexes, Gs-9667 had the lowest mean value of RMSD. Additionally, RMSF analysis revealed the high fluctuations at the looped regions within the nsp16 structure, particularly around the amino acids 28, 140, and 260. Even though the average fluctuations of the eight protein-drug complexes were very similar, the 2 0 -O-MTase-Hesperidin complex showed the most significant fluctuations. Furthermore, the conserved K-D-K-E (K46, D130, K170, and E203) motif remained remarkably stable throughout the MD simulation. Our MD simulation analysis also found an average of five hydrogen bonds present between nsp16 protein and Hesperidin. Structural analysis indicated that Hesperidin might form hydrogen bonds with K46 and K170 of the KDKE motif, as well as N43 and L100, which are key residues for SAM binding in nsp16. Recent virtual screening studies also found that Hesperidin may disrupt the interaction of SARS-CoV-2 Spike protein and its host cellular receptor ACE2 (Angiotensin-converting enzyme 2) (Haggag et al., 2020; . It will be interesting to determine how effectively Hesperidin binds nsp16 of SARS-CoV-2, and to demonstrate its antiviral activity in cellular and animal models related to COVID-19.

In summary, we have analyzed structural features and performed functional determination of 2 0 -O-MTase nsp16 of SARS-CoV-2, which causes the devastating COVID-19. The nsp16 has a druggable pocket, which is associated with substrate binding and enzymatic activity. Using high-throughput virtual screening, molecular docking, and MD simulation approaches, we identified 17 therapeutic drugs, including several clinically approved drugs that have the potential to inhibit 2 0 -O-MTase of SARS-COV-2. The identified top candidates, including Hesperidin, Rimegepant, Gs-9667, and Sonedenoson, merit further testing and repurposing experimentally and/or clinically for COVID-19 treatment.

The primary amino acid sequences of seven human CoV nsp16 proteins and human CMTR1 were retrieved from the NCBI (National Center for Biotechnology Information) Database, and Supplementary Table 1 lists the detailed genome and protein ID for these 2 0 -O-MTases. The protein sequences were aligned using Clustal Omega (https://www. ebi.ac.uk/Tools/msa/clustalo/), and was presented with the ESPript 3.0 program (http://espript.ibcp.fr/ESPript/ESPript/). Crystal structures of three CoV nsp16s and human CMTR1 were obtained from RCSB Protein Data Bank (Berman et al., 2000) , and analyzed with PyMOL, UCSF Chimera and Protein-Ligand Interaction Profiler programs (Pettersen et al., 2004; Salentin et al., 2015) . Structural alignment and comparison were performed using TM-align algorithm and Dali server (Holm, 2020; Zhang & Skolnick, 2005) .

Crystal structure (PDB: 6W4H) of SARS-CoV-2 nsp16 with 1.8 Å resolution was prepared with UCSF Chimera programs for virtual screening against a Drugs-lib that contains 7,173 stereoisomers corresponding to 4,574 "approved" drugs in MTiOpenScreen service (Lagarde et al., 2018 (Lagarde et al., , 2019 . The Drugs-lib is generated from four compound databases of the "drug" subset of the ChEMBL database, the "approved" subset of DrugBank, the DrugCentral database, and the "approved" SuperDrug2 database (Lagarde et al., 2019) . For the MTiOpenScreen Vina docking, the (x, y, z) grid center coordinates used for SAM binding pocket of 6W4H are (83.7, 14.8, 27.7) , and the size of the search space was set to 20 Å x 20 Å x 20 Å. The MTiOpenScreen screening was repeated three times, and the compounds shared across all three runs were analyzed by DataWarrior, PyMOL, and UCSF Chimera programs. Fragfp descriptors in the DataWarrior were used to compute the similarity between each pair of compounds (Sander et al., 2015) . The targets and pathways of identified candidate compounds were analyzed by querying Probes and Drugs Portal (Skuta et al., 2017) . The interaction of candidate drugs with key residues for nsp16 2 0 -O-MTase function were analyzed by PyMOL, LIGPLOT, and PLIP (protein-ligand interaction profiler) (Salentin et al., 2015; Wallace et al., 1995) .

Ligand flexible docking simulation was performed using the RosettaLigand program that employs the Monte Carlo minimization protocol (Davis & Baker, 2009; Lemmon & Meiler, 2012) . In the first, low-resolution stage, the ligand was placed near the SAM binding pocket of nsp16, where it was perturbed 50 times and rotated 1,000 times in a random direction. The best-scoring models were filtered by rootmean-square deviation to eliminate near-duplicates and one of the remaining models was selected at random. To ensure that the ligand did not "walk away" from the protein, the ligand was moved towards the protein. The second, high-resolution stage employed the Monte Carlo minimization protocol in which the ligand position and orientation were randomly perturbed by a small deviation (0.1 Å and 3 ); receptor side chains were repacked using a rotamer library; the ligand position, orientation, and torsions and protein side-chain torsions were simultaneously optimized using quasi-Newton minimization; and the end result was accepted or rejected based on the Metropolis criterion. Scoring used the full-atom Rosetta energy function with softened van der Waals repulsion (Davis & Baker, 2009; Lemmon & Meiler, 2012) . The full repack made 1,000 random rotamer substitutions at random positions and accepted or rejected each on the Metropolis criterion. Rotamer trials chose the single best rotamer at a random position in the context of the current state of the rest of the system, with the positions visited once each in random order. Six cycles of rotamer trials and a full repack after every three cycles were performed. The third and final stage was a more stringent, gradient-based minimization of the ligand position, orientation, and torsions and receptor torsions for both side chains and the backbone. Scoring used the same Rosetta energy function, but with a hard-repulsive van der Waals potential, which created a more rugged energy landscape that was better at discriminating native from non-native binding modes (Davis & Baker, 2009; Lemmon & Meiler, 2012) . The entire structure was then scored using the partial covalent interactions energy function, which was developed to score H-bonds more accurately.

To make the docking results mimic the physiological state of protein molecules, we performed MD simulations on PDB 6W4H, which contained both nsp16 and nsp10, with the eight selected inhibitors and natural ligand SAM. The MD simulations were executed by the GROningen MAchine for Chemical Simulations (GROMACS) version 2020.1 with GROMOS96 43a1 force field parameters (Chiu et al., 2009) . The topology of the drug molecules was created using PRODRG webserver (Schuttelkopf & van Aalten, 2004) . The simulation of all the proteins and protein-drug complexes were run for a period of 100 ns. To make the system electrostatically neutral, counter ions were added to the proteindrug complexes. The complexes were solvated within a 10 nm SPC/E (extended simple point charge) water cube . The protein-drug complexes were minimized in multiple steps using the steepest descent method, where minimizations of the whole system, water cube, and non-heavy atoms of the complexes were accomplished. The entire systems were then progressively heated up to 300 K on a time scale of 100 ps. The equilibration steps were performed in two different phases, one with constant pressure and temperature (NPT) and the other with steady volume and temperature (NVT). Various structural parameters, like Root Mean Square Deviation (RMSD), Root Mean Square Fluctuations (RMSF), Radius of Gyration (Rg), Principal Component Analysis (PCA) based on essential dynamics (ED) approach, Inter-molecular Hydrogen Bonding (H-bonding), and Solvent-Accessible Surface Area (SASA) were calculated as a function of time to explore the structural behavior of the proteins and protein-drug complexes.

The binding affinity between the protein and drug molecules were calculated using molecular mechanics/Poisson-Boltzmann surface area (MM-PBSA) method [84] . MM-PBSA utilizes the effects of thermal averaging with a force field/ continuum solvent model to post-process a series of representative snapshots from MD trajectories (Poli et al., 2020) . The method expresses the binding affinity (DG bind) as the difference between the free binding energy of the complex and the free binding energy of the receptor plus the ligand. In this study, the binding free energy was calculated as the Figure 7 . The per-residue free energy contribution spectrums of nsp16 in four compound-protein complexes. Only residues contributing above a ± 10 kcal/mol threshold were colored and labeled.

average of the MM-PBSA obtained from 400 equally distributed snapshots of MD trajectory of each protein-ligand complex after removing the water molecules by using GMXPBSA 2.1 tool available in GROMACS (Paissoni et al., 2015) . A g_mmpbsa tool was used calculate the contribution of each residue to the total binding free energy (Kumari et al., 2014) .

No potential conflict of interest was reported by the authors.

This work was partially supported by grants from the Department of Defense (DoD) Breast Cancer Program BC161536, Elsa U. Pardee Foundation, DMC Foundation and Molecular Therapeutics Program of Karmanos Cancer Institute to Dr. Zeng-Quan Yang; and by funding from Susan G. Komen GTDR14299438 and Wayne State University Graduate School Dean Mathur Fellowship to Morenci Manning.

Molecular mechanisms of action of hesperidin in cancer: Recent trends and advancements

Binding of the methyl donor S-adenosyl-l-methionine to Middle East respiratory syndrome coronavirus 2'-O-methyltransferase nsp16 promotes recruitment of the allosteric activator nsp10

RNA ribose methylation (2'-O-methylation): Occurrence, biosynthesis and biological functions

Antidiabetic potential of flavonoids from traditional Chinese medicine: A review

Characterization of hMTr1, a human Cap1 2'-O-ribose methyltransferase

Calcitonin gene-related peptide receptor antagonists: New therapeutic agents for migraine

The protein data bank

A new insight into protein-protein interactions and the effect of conformational alterations in PCNA

Identification of bioactive molecules from tea plant as SARS-CoV-2 main protease inhibitors

Modomics: A database of RNA modification pathways

In vitro reconstitution of SARS-coronavirus mRNA cap methylation

RNA methylation under heat shock control

RNA methyltransferases involved in 5' cap biosynthesis

Structure-based virtual screening and experimental validation of the discovery of inhibitors targeted towards the human coronavirus nucleocapsid protein

Biochemical and structural insights into the mechanisms of SARS coronavirus RNA ribose 2'-O-methylation by nsp16/nsp10 protein complex

An improved united atom force field for simulation of mixed lipid bilayers

The species Severe acute respiratory syndrome-related coronavirus: Classifying 2019-nCoV and naming it SARS-CoV-2

RosettaLigand docking with full ligand and receptor flexibility

Crystal structure and functional analysis of the SARS-coronavirus RNA cap 2'-O-methyltransferase nsp10/ nsp16 complex

Hesperidin attenuates influenza A virus (H1N1) induced lung injury in rats through its anti-inflammatory effect

TRK inhibitors in TRK fusion-positive cancers

Potential drugs targeting early innate immune evasion of SARS-coronavirus 2 via 2'-O-methylation of viral RNA

Mutations in the FTSJ1 gene coding for a novel S-adenosylmethionine-binding protein cause nonsyndromic X-linked mental retardation

Emerging adenosine receptor agonists: An update

Gene expression profiling of the cellular transcriptional network regulated by alpha/beta interferon and its partial attenuation by the hepatitis C virus nonstructural 5A protein

Dynamic docking: A paradigm shift in computational drug discovery

Making the mark: The role of adenosine modifications in the life cycle of RNA viruses

Defects in tRNA anticodon loop 2'-O-methylation are implicated in nonsyndromic Xlinked intellectual disability due to mutations in FTSJ1

Active site in RrmJ, a heat shock-induced methyltransferase

Is hesperidin essential for prophylaxis and treatment of COVID-19 Infection? Medical Hypotheses

Knowledge-based structural models of SARS-CoV-2 proteins and their complexes with potential drugs

Using Dali for protein structure comparison

Innate immune restriction and antagonism of viral RNA lacking 2‫-׳‬O methylation

Regulation and function of CMTR1-dependent mRNA cap methylation

Historical and current adenosine receptor agonists in preclinical and clinical development

Epitranscriptomic code and its alterations in human disease

Targeting SARS-CoV-2: A systematic drug repurposing approach to identify promising inhibitors against 3C-like proteinase and 2'-O-ribose methyltransferase

Docking and scoring in virtual screening for drug discovery: Methods and applications

g_mmpbsa-a GROMACS tool for high-throughput MM-PBSA calculations

A free web-based protocol to assist structure-based virtual screening experiments

Online structure-based screening of purchasable approved drugs and natural compounds: Retrospective examples of drug repositioning on cancer targets

Discovery of a selective catalytic p300/CBP inhibitor that targets lineage-specific tumours

Assignment of 2'-O-methyltransferases to modification sites on the mammalian mitochondrial large subunit 16 S ribosomal RNA (rRNA)

Mitochondrial ribosomal RNA (rRNA) methyltransferase family members are positioned to modify nascent rRNA in foci near the mitochondrial DNA nucleoid

Rosetta Ligand docking with flexible XML protocols

Genome-wide CRISPR screen identifies host dependency factors for influenza A virus infection

The pivotal regulatory landscape of RNA modifications

Rimegepant, an oral calcitonin gene-related peptide receptor antagonist, for migraine

Intellectual disability-associated gene ftsj1 is responsible for 2'-O-methylation of specific tRNAs

Structure-based virtual screening: From classical to artificial intelligence

Pan-cancer analysis of RNA methyltransferases identifies FTSJ3 as a potential regulator of breast cancer progression

A first in man, dose-finding study of the mTORC1/mTORC2 inhibitor OSI-027 in patients with advanced solid malignancies

Coronavirus nonstructural protein 16: Evasion, attenuation, and possible treatments

Attenuation and restoration of severe acute respiratory syndrome coronavirus mutant lacking 2'-o-methyltransferase activity

Epitranscriptomic marks: Emerging modulators of RNA virus gene expression

GMXPBSA 2.1: A GROMACS tool to perform MM/PBSA and computational alanine scanning

Structure-based drug designing and immunoinformatics approach for SARS-CoV-2

Coronaviruses post-SARS: Update on replication and pathogenesis

UCSF Chimera-a visualization system for exploratory research and analysis

Molecular docking: Shifting paradigms in drug discovery

Application of MM-PBSA methods in virtual screening

Pathological role of a point mutation (T315I) in BCR-ABL1 protein-A computational insight

FTSJ3 is an RNA 2'-O-methyltransferase recruited by HIV to avoid innate immune sensing

Dynamic RNA modifications in gene expression regulation

PLIP: Fully automated protein-ligand interaction profiler

Pharmacological properties of MK-3207, a potent and orally active calcitonin gene-related peptide receptor antagonist

DataWarrior: An open-source program for chemistry aware data visualization and analysis

Entrectinib for the treatment of metastatic NSCLC: Safety and efficacy

PharmaGist: A webserver for ligand-based pharmacophore detection

PRODRG: A tool for highthroughput crystallography of protein-ligand complexes

Rimegepant: First Approval

Insights into RNA synthesis, capping, and proofreading mechanisms of SARScoronavirus

Where, when, and how: Context-dependent functions of RNA methylation writers, readers, and erasers

Probes &Drugs portal: An interactive, open data resource for chemical biology

The compromise of virtual screening and its impact on drug discovery

Protein structure-based drug design: From docking to molecular dynamics

Structural analysis of human 2'-O-ribose methyltransferases involved in mRNA cap structure formation

Reduction of free fatty acids, safety, and pharmacokinetics of oral GS-9667, an A(1) adenosine receptor partial agonist

Genomic analysis of the host response to hepatitis C virus infection

SARS-CoV ORF1b-encoded nonstructural proteins 12-16: Replicative enzymes as antiviral targets

Targeting SARS-COV-2 non-structural protein 16: A virtual drug repurposing study

Potential anti-inflammatory effects of hesperidin from the genus citrus

AutoDock Vina: Improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading

GROMACS: Fast, flexible, and free

Identification of promising drug candidates against NSP16 of SARS-CoV-2 through computational drug repurposing study

Structural basis of RNA cap modification by SARS-CoV-2

DoGSiteScorer: A web server for automatic binding site prediction, analysis and druggability assessment

LIGPLOT: A program to generate schematic diagrams of protein-ligand interactions

Analysis of therapeutic targets for SARS-CoV-2 and discovery of potential drugs by computational methods

A new coronavirus associated with human respiratory disease in China

TM-align: A protein structure alignment algorithm based on the TM-score

A novel coronavirus from patients with pneumonia in China

Ribose 2'-O-methylation provides a molecular signature for the distinction of self and non-self mRNA dependent on the RNA sensor Mda5