Computational guided drug repurposing for targeting 2'-O-ribose methyltransferase of SARS-CoV-2.

AIMS: The recent outbreak of pandemic severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2) has led the world towards a global health emergency. Currently, no proper medicine or effective treatment strategies are available; therefore, repurposing of FDA approved drugs may play an important role in overcoming the situation.
MATERIALS AND METHODS: The SARS-CoV-2 genome encodes for 2-O-methyltransferase (2'OMTase), which plays a key role in methylation of viral RNA for evading host immune system. In the present study, the protein sequence of 2'OMTase of SARS-CoV-2 was analyzed, and its structure was modeled by a comparative modeling approach and validated. The library of 3000 drugs was screened against the active site of 2'OMTase followed by re-docking analysis. The apo and ligand-bound 2'OMTase were further validated and analyzed by using molecular dynamics simulation. KEY
FINDINGS: The modeled structure displayed the conserved characteristic fold of class I MTase family. The quality assessment analysis by SAVES server reveals that the modeled structure follows protein folding rules and of excellent quality. The docking analysis displayed that the active site of 2'OMTase accommodates an array of drugs, which includes alkaloids, antivirals, cardiac glycosides, anticancer, steroids, and other drugs. The redocking and MD simulation analysis of the best 5 FDA approved drugs reveals that these drugs form a stable conformation with the 2'OMTase. SIGNIFICANCE: The results suggested that these drugs may be used as potential inhibitors for 2'OMTase for combating the SARS-CoV-2 infection.

Introduction

The recent outbreak of severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2) in Wuhan, China, in December 2019, has infected over 69,16,233 and death of over 4,00,135 (Dated June, 7, 2020) worldwide since its outbreak. The number of cases is increasing every hour, and due to this reason, the World Health Organization (WHO) has declared it as a global public health emergency.

The SARS-CoV-2 is a positive-stranded RNA virus with a genome size of 29.9 kb encoding 12 open reading frames (ORFs) and two untranslated regions of 254 and 229 bp at 5′ and 3′ ends, respectively [1]. The molecular organization of the coronavirus genome in 5′ → 3′ direction encodes several proteins. The genome includes non-structural and structural proteins. The structural proteins include spike (S), envelope (E), matrix (M), and nucleocapsid (N) proteins [1]. The ~66% of the genome at the 5′end contains two overlapping ORFs (ORF1a and ORF1b) encoding polyproteins 1a (PP1a) and 1ab (PP1ab). These polyproteins are processed into 16 functional non-structural proteins by papain-like protease and 3C-like protease, which form the replication transcription complex (RTC). These non-structural proteins include proteases and RNA processing enzymes [[2], [3], [4]]. These RNA processing enzymes are uncommon and include exoribonuclease, guanine N7-methyltransferase (N7-MTase), and 2′O-methyltransferase. Approximately, 33% of the genome at 3′ end encodes four major structural proteins S, E, M, and N. These structural proteins play a crucial role in the formation of the new virus particle and its release [5].

Most viral and eukaryotic cellular mRNA containing polyadenine region at 3′ terminal and a 5′terminal cap structure of N7-methyl-guanine moiety connected through 5′-5′ triphosphate bridge to the first transcribed nucleotide [6,7]. The 5′ terminal cap structure of mRNA plays a key role in protecting its degradation by 5′-exonuclease and efficient translation. The capping of mRNA (Cap0) is catalyzed by RNA triphosphatase, a guanylyltransferase and guanine-N7-methyltransferase (N7-MTase) in a sequential reaction. The capping of most viruses and eukaryotic mRNA is further methylated at the 2-OH position of the first nucleotide by ribose 2′-O-methyltransferase (2′MTase) enzyme to form the cap1 structure. The 2′O-methylation of viral mRNA cap also helps in protecting RNA recognition as “nonself” by host innate immune system, which makes it a potential target for antiviral drug development [7].

The design of antiviral drugs with proven efficacy or vaccines against SARS-CoV-2 is limited due to the poor understanding of its molecular mechanisms of infection. Therefore, there is a need to design an inhibitor for blocking the SARS-CoV-2 replication within the cell. The repurposing of the approved drugs is an alternate way to find out a solution. Recently, clinicians have used the antimalarial and anti-HIV drugs for the treatment of SARS-CoV-2; however, the possible mechanism involved in its inhibition is not known [[8], [9], [10]].

In the present study, the 2′-O-methyltransferase (2′OMTase) structure was modeled using a comparative modeling approach and further validated. Besides, the FDA approved drugs were screened against 2′OMTase for selecting the potential candidate for drug repurposing using virtual drug screening followed by redocking analysis using Autodock tool. The selected best FDA approved drugs were further validated by performing the molecular dynamics simulation.

Material and methods

Amino acid sequence retrieval and analysis

The protein sequence encoding ORF1ab polyprotein of Indian isolates of SARS-CoV-2 (GenBank Accession ID: QIA98582.1) constituting 7096 amino acid residues was retrieved from the NCBI Database (https://www.ncbi.nlm.nih.gov/protein/1809484466). The retrieved sequence was pairwise aligned with the 2′-O-ribose methyltransferase (NCBI RefSeq ID: YP_009725311.1) encoding region of SARS-CoV-2 isolate from Wuhan to identify the 2′OMTase. PSI-BLAST webserver was used to search for 2′OMTase domains from the different strains across the SwissProt database. 2′OMTase sequence was aligned with its homologous characterized domains with a query coverage of 100% displaying the sequence identity in the range of 53%–93% by T-Coffee Multiple sequence alignment web server (http://tcoffee.crg.cat/apps/tcoffee/do:expresso). Furthermore, the results were visualized and analyzed for the identification of conserved and semi-conserved amino acid residues by Esprit 3.0 webserver (http://espript.ibcp.fr/ESPript/ESPript/).

Comparative modeling, refinement and structure assessment of 2′OMTase

The structure modeling of 2′OMTase was performed by using the comparative modeling approach through Modeller v9.23 Program [11]. The homologous protein structures of 2′OMTase with PDB ID: 3R24 [12,13], 2XYR [12,13] from human SARS-CoV and 2′OMTase with PDB ID 5YN5 from Human betacoronavirus 2c EMC/2012 [Unpublished] were retrieved from Protein Data Bank (PDB). The structures mentioned above were aligned by running ‘salign( )’ command to generate the structure-based multiple sequence alignment. After that, the generated Multiple Sequence Alignment (MSA) was further aligned with the 2′OMTase query sequence to cover the whole sequence for model building. A total of 20 independent models of 2′OMTase were built by running the ‘automodel’ command, and the best two models displaying least DOPE scores were selected for energy minimization by YASARA webserver. The energy minimized models of 2′OMTase were evaluated by using several programs viz, ProCheck [14], ERRAT [15] and VERIFY3D [16] tools available at SAVES server (https://services.mbi.ucla.edu/SAVES/) and RAMPAGE (http://mordred.bioc.cam.ac.uk/~rapper/rampage.php) server. The ProSA web server (https://prosa.services.came.sbg.ac.at/prosa.php) was used to estimate the statistical Z-score deviation of the modeled structure from the high-resolution structures deposited in PDB. The best-selected energy minimized 2′OMTase modeled structure was used for further analysis. The topology diagram showing the secondary structures present in 2′OMTase was generated by using PDBSum (http://www.ebi.ac.uk/thornton-srv/databases/pdbsum/Generate.html).

Screening of FDA drugs approved against 2′OMTase

To identify the potential antivirals drugs candidate for repurposing against 2′OMTase, the PyRx 0.8 virtual screening tool was used. The library includes 3000 FDA approved drugs containing antivirals, antibiotics, protease inhibitors, anticancer, membrane transport and ion channel inhibitors, DNA damage and repair inhibitors etc. All the FDA approved drugs were downloaded from the zinc database and DrugBank and converted into 3D format for high-throughput virtual screening. All the FDA approved drugs were energy minimized and converted to autodock vina supporting format by adding quaternion and torsions by following the PyRx 0.8 protocol. The prepared FDA approved drugs were screened against the catalytic site of SARS-CoV-2′OMTase protein using AutoDock Vina program compiled in PyRx 0.8 by setting the grid center points at X = −12.32, Y = 0.34, Z = 9.52 and box dimensions were set as 32.15 Å × 32.34 Å × 35.51 Å with the exhaustiveness of 8. The drug molecules showing the free energy of binding (ΔG) between −10 kcal/mol −8.7 kcal/mol were selected for further analysis. The compound sinefungin, binds in the S-adenosylmethionine binding and resulted in the regulation of RNA capping enzymes of SARS-CoV, therefore, it was used as reference [13]. The free energy of binding of −8.3 kcal/mol for sinefungin was taken as reference.

Molecular docking analysis of 2′OMTase with selected FDA approved drugs

The selected best FDA approved drugs obtained from virtual screening were redocked using AutoDock v4.2.6 software integrated with MGL Tools 1.5.6 (http://mgltools.scripps.edu/). The ligands and 2′OMTase were prepared and converted into the PDBQT format. The Grid box of 70x70x70 Å size in x, y and z-direction with default grid point spacing was arranged around the active site by using the grid centering option. The conformational search and docking simulation for different ligands were executed using 2000 cycles of the Lamarckian Genetic Algorithm (LGA). The rest of the parameters were set as default. These 2000 conformations of ligands were ranked and clustered based on the free energy of binding. The best conformation was selected based on the free energy of binding for protein-ligand interaction analysis and estimation of inhibitory constant. Top five FDA approved drugs were further used for analysis. The protein-ligand complexes were further analyzed by PyMOL for the identification of hydrophobic contacts and polar interactions. The 2D interaction of protein and ligand was generated by using the LigPlot+ program (http://www.ebi.ac.uk/thornton-srv/databases/pdbsum/Generate.html).

Molecular dynamic simulation of the apo and ligand-bound form of 2′OMTase

The conformational stability and behavioral dynamics of apo and ligand-bound 2′OMTase structures were determined by molecular dynamic (MD) simulation. GROMOS96 43a1 force field compiled with GROMACS v5.1.4 software program installed in the High-Performance Computing facility (Param-Ishan) available at Indian Institute of Technology Guwahati, India was used for MD simulation. The topology of FDA approved drug viz. Sinefungin, dihydroergotamine, digitoxin, irinotecan, and teniposide were generated by using PRODRG server (http://prodrg1.dyndns.org/submit.html). The apo form and ligand-bound complex was placed in a triclinic box with a distance of 1.2 nm and solvated with a single point charge (SPC) waters. The system was neutralized by adding 2 Cl− ions as counterions, and the system was energy minimized by using the Steepest descent method using cut-off up to 239 kcal/mol. The system was equilibrated in two steps NVT ensemble (constant number of particles, volume and temperature) followed by NPT ensemble (constant number of particles, pressure and temperature) for 500 ps each. The equilibrated system was simulated for 50 ns to run the final MD simulation, and the trajectories were recorded at 10 ps interval. The final MD run was analyzed by using gmx rms for the estimation of root-mean-square deviation (RMSD), gmx rmsf for room mean square function (RMSF) calculation. The radius of gyration (Rg) and solvent accessible surface area (SASA) were determined by using gmx gyrate, and gmx sasa commands, respectively. The nonbonded interaction energy between protein and ligand was estimated using energygrps functions compiles in gromacs v5.14. It helps in estimating the reliability of non-bonded interactions. The XMGRACE software (http://plasma-gate.weizmann.ac.il/Grace/) was used for data plotting and analysis.

Results

Sequence analysis of 2′OMTase

The pairwise sequence alignment analysis of ORF1ab polyprotein of Indian isolate of SARS-CoV-2 with isolate from Wuhan revealed that both the sequences of the 2′-O-ribose methyltransferase domain are identical without any mutation. The PSI-BLAST analysis of 2′OMTase displayed the amino acid sequence identity with previously characterized 2′OMTase of CoV isolated from different strains (Table 1 ). The 2′OMTase sequences from different strains mentioned in Table 1 were retrieved from Uniprot Database and aligned using the T-Coffee program revealed that the catalytic amino acid residues Lys46, Asp130, Lys170, and Glu203 are conserved in all the strains (Fig. 1 ).

Table 1

PSI-BLAST analysis of SARS-CoV 2′OMTase with homologous structural enzymes obtained from SwissProt database.

UniProt ID	Virus	Query cover (%)	% identity
P0C6X7	Human SARS coronavirus	100.00	93.30
K9N7C7	Middle East respiratory syndrome virus	100.00	66.30
P0C6W3	Bat coronavirus HKU4	100.00	66.60
P0C6W5	Bat coronavirus HKU9	100.00	65.10
P0C6X9	Murine coronavirus (strain A59)	100.00	65.20
P0C6X5	Human coronavirus NL63	100.00	58.40
P0C6X1	Human coronavirus 229E	100.00	57.00
P0C6Y5	Porcine transmissible gastroenteritis coronavirus	100.00	55.40
Q98VG9	Feline coronavirus (FIPV WSU-79/1146)	100.00	55.40
P0C6Y4	Porcine epidemic diarrhea virus (CV777)	100.00	57.30
P0C6Y1	Avian infectious bronchitis virus (strain Beaudette)	100.00	53.00
P0C6Y3	Avian infectious bronchitis virus (strain M41)	100.00	52.70

Fig. 1

Multiple Sequence Alignment (MSA) of 2′OMTase with characterized 2′OMTase from different strains of SARS-CoV, as given in Table 1. The conserved amino acid residues are displayed in the red color background and semi-conserved residues are shown in red color. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Structure modeling of 2′OMTase and validation

The comparative modeled structure of 2′OMTase (Fig. 2A) displayed the characteristic fold comprising seven stranded β-sheets surrounded by α-helices and loops (Fig. 2A). A similar type of structural organization is conserved and reported for class I MTase family [12,13]. The different programs available on the SAVES server were used to assess the quality of 2′OMTase energy minimized modeled structure. The Ramachandran plot analysis by Procheck server displayed that 266 amino acid residues i.e., non-glycine and proline, follow the stereochemical properties. The analysis showed that 91.4% amino acid residues of 266 lies in the favored region, and 8.6% amino acid residues accommodate in the allowed region (Fig. 2B). The absence of amino acid residues in the generously allowed and disallowed regions of the Ramachandran plot demonstrated that the modeled structure of 2′OMTase follows the stereochemical property rule and follows all the possible dihedral, phi (ϕ) and psi (ψ) angle values. The modeled structure on further evaluation by the VERIFY3D plot, revealed that the 92.28% of amino acid residues have an overall average 3D-1D score ≥ 0.2 (Fig. 2C). Similarly, the 2′OMTase energy minimized structure assessed by ERRAT Plot displayed the quality factor of 96.20%, which further confirmed that the modeling of 2′OMTase structure is error-free (Fig. 2D). Similarly, ProSA analysis of 2′OMTase displayed the Z score of −7.4, confirming the modeled structure accommodated in the X-ray zone (Fig. 2E). PDBSum analysis of the modeled structure of 2′OMTase displayed that it contains 13 α-helices and 13 β-strands as secondary structure elements (Fig. 2F). The quality assessment parameters of the modeled structure of 2′OMTase were superior to the experimentally determined homologous crystal structure of 2′OMTase, which reveals that the selected modeled structure is the best model for further studies. The 2′OMTase and its homologous proteins were subjected to CastP analysis for the determination of active site volume and to predict the open or close conformation. The active site volume of 2′OMTase from SARS-CoV-2 was found to be 526.70 Å3. The active site volume of homologous 2′OMTase for PDB ID 3R24_A was 946.75 Å3 and for PDB ID 2XYR was 537.38 Å3. The Active site volume of 2′OMTase from SARS-CoV-2 matched well with the homologous protein of SARS-CoV and confirmed that the modeled protein has closed conformation.

Fig. 2

Three-dimensional structure modeling of 2′OMTase and validation analysis. A) 3D structure B) Ramachandran Plot analysis, C) VERIFY 3D showing the threshold score >0.2, D) ERRAT Plot analysis showing the quality of the 2′OMTase modeled structure E) ProSA plot showing Z-score and F) 2′OMTase Topology diagram.

Screening of FDA approved drugs against 2′OMTase

The high throughput virtual screening analysis of FDA approved drugs library against the 2′-O-methyltransferase (2′OMTase) helps in identifying and selecting the potential FDA approved drugs as an inhibitor or antiviral compounds. Top twenty two molecules displaying the free energy of binding in the range of −10.00 to −8.7 kcal/mol were selected and considered as top hit compounds and are reported in Table 2 . Each ligand screened through virtual screening by PyRx displayed nine different conformations, and the ligands forming interaction with catalytic tetrad were selected for further analysis. These selected top hits included antiviral, alkaloids based drugs, cardiac glycoside, anticancer drugs, steroid-based drugs, and other drugs, etc. These selected FDA approved drugs were further employed for re-docking analysis with the 2′OMTase to discover the potential inhibitor(s) by using Autodock tool.

Table 2

Best FDA approved drugs obtained from Virtual screening.

Ligands	Binding energy(kcal/mol)
Alkaloids based drugs
Ergotamine	−10
Dihydroergotamine	−9.5
Antiviral
Saquinavir	−8.7
Indinavir	−8.7
Cardiac glycoside
Digitoxin	−9.4
Deslanoside	−8.9
Chlorthalidone	−8.7
Anticancer compounds
Dactinomycin	−9.9
Irinotecan	−9.4
Nilotinib	−9.4
Plicamycin	−9.4
Methotrexate	−8.8
Eribulin	−8.7
Lanreotide	−8.7
Steroid based drugs
Fluocinolone acetonide	−8.9
Desoximetasone	−8.8
Triamcinolone acetonide	−8.7
Betamethasone	−8.7
Other drugs
Tadalafil	−9
Zafirlukast	−8.9
Paliperidone	−8.8
Ertapenem	−8.7

Re-docking analysis of 2′OMTase with identified FDA approved drugs

The re-docking analysis of the top twenty selected screened FDA approved drugs was carried out by AutoDock 4.2.6. The molecular docking analysis revealed that the top five FDA approved drugs viz. sinefungin, dihydroergotamine, digitoxin, irinotecan, and teniposide displayed the free energy of binding in the range of −9.8 kcal/mol to −7.8 kcal/mol and predicted inhibitory constant in the range of 20 nM to 3120 nM (Table 3 ). The amino acid residues involved in the formation of hydrogen bonds, other interactions, the free energy of binding and inhibitory constant are shown in Table 3. The protein-ligand complex of selected FDA approved drugs and 2′OMTase were further analyzed by using PyMOL v2.3 software and LigPlot+ software.

Table 3

FDA approved drugs obtained from re-docking by Autodock.

Ligands	Binding energy(kcal/mol)	Predicted inhibitory constant (nM)	Amino acid involved in H-bond/polar contacts (with residue no.)	Amino acid residues within 4 Å of active site and involved in hydrophobic/Van der Waals interaction
Control
Sinefungin	−7.8	3120.05	Gly71, Asp130, Tyr132, Lys170	Tyr47, Gly73, Ser74, Asn99, Leu100, Met131, Asp133, Pro134, Phe149
Alkaloids based drugs
Dihydroergotamine	−9.3	279.49	Tyr47, Gly71, Lys170	Asn43, Gly73, Ser74, Asp75, Lys76, Pro80, Asp99, Leu100, Asn101, Asp130, Met131, Tyr132, Pro134, Glu203
Cardiac glycoside
Digitoxin	−9.1	162.91	Lys24, Tyr30, Tyr132 His174	Asn43, Gly71, Gly73, Ser74, Asp75, Asp99, Leu100, Asp130, Met131, Asp133, Pro134, Thr136, Phe149, Lys170, Glu173, Ser202, Glu203
Anticancer compounds
Irinotecan	−9.3	20.3	Tyr47, Ser200	Tyr30, Gly31, Asp32, Ser33, Met42, Asn43, Lys46, Gly71, Gly73, Ser74, Asp99, Leu100, Asp130, Met131, Tyr132, Lys170, Val197, Asn198, Ser201
Teniposide	−8.8	450.3	Asp99, Leu100, Asp130, Lys137, Glu203	Gly71, Gly73, Ser74, Asn101, Met131, Tyr132, Pro134, Lys170, Thr172

The drug, sinefungin, a natural nucleoside, interacts with 2′OMTase through four hydrogen bonds formed by Gly71, Asp130, Tyr132, and Lys170 (Fig. 3A). The sinefungin interacts with Tyr47, Gly73, Ser74, Asn99, Leu100, Met131, Asp133, Pro134, and Phe149 amino acid residues and forms hydrophobic/Van der Waals interaction (Fig. 3A). The amino acid residue Asp130 and Lys170 forming hydrogen bonds with sinefungin, are the part of catalytic tetrad and plays a key role in the formation of the catalytic cleft.

Fig. 3

Molecular docking analysis of 2′OMTase with best FDA approved drugs. A) Sinefungin B) Digitoxin, C) Dihydroergotamine, D) Irinotecan and E) Teniposide.

Similarly, digitoxin, a cardiac glycoside interacts with 2′OMTase by forming four hydrogen bonds by Lys24, Tyr30, Tyr132 and His174 followed by the formation of hydrophobic interaction through Asn43, Gly71, Gly73, Ser74, Asp75, Asp99, Leu100, Asp130, Met131, Asp133, Pro134, Thr136, Phe149, Lys170, Glu173, Ser202 and Glu203 (Fig. 3B). The amino acid residues of catalytic tetrad viz., Asp130, Lys170, and Glu203 involved in hydrophobic interactions revealed that these residues are further involved in stabilizing the ligand in the active site region (Fig. 3B).

The alkaloid drug, dihydroergotamine, forms hydrogen bonds with Tyr47, Gly71, and Lys170 and forms hydrophobic interaction with Asn43, Gly73, Ser74, Asp75, Lys76, Pro80, Asp99, Leu100, Asn101, Asp130, Met131, Tyr132, Pro134 and Glu203 of 2′OMTase (Fig. 3C). The three amino acid residues of catalytic tetrad stabilize the drug in the catalytic cleft by hydrogen bond and hydrophobic interactions.

Anticancer drug irinotecan (Fig. 3D) forms hydrogen bonds with Tyr47 and Ser200 and involved in hydrophobic interactions through Tyr30, Gly31, Asp32, Ser33, Met42, Asn43, Lys46, Gly71, Gly73, Ser74, Asp99, Leu100, Asp130, Met131, Tyr132, Lys170, Val197, Asn198 and Ser201 of 2′OMTase. Teniposide interacts with Asp99, Leu100, Asp130, Lys137, and Glu203 and forms hydrogen bond (Fig. 3E). Besides, the amino acid residues Gly71, Gly73, Ser74, Asn101, Met131, Tyr132, Pro134, Lys170, and Thr172 form hydrophobic interaction. These drugs are stabilized in the catalytic cleft by catalytic residues Asp130, Lys170, and Glu203.

The selected FDA approved drugs completely accommodated the catalytic cleft and interaction stabilized mainly by three amino acid residues viz. Asp130, Lys170, and Glu203 of the catalytic tetrad. This may result in preventing the 5′ Cap modification by inhibiting the methyltransferase enzyme activity.

The molecular dynamics (MD) simulation analysis of both apo and ligand complex was analyzed for the 50 ns to understand the stability and conformational dynamic behavior. The estimated root mean square deviation (RMSD) analysis of apo 2′OMTase (Black color) displayed fluctuation up to eight ns, and thereafter it became stable at 0.24 nm till 50 ns (Fig. 4A). The drug sinefungin bound 2′OMTase displayed the RMSD at 0.30 nm till 25 ns, but after that, the sudden increase was observed until 50 ns. Similarly, digitoxin bound 2′OMTase displayed the initial RMSD at 0.32 nm, and after that, a small increase from 0.32 nm to 0.36 nm was observed up to 50 ns. The drug dihydroergotamine bound 2′OMTase showed an average RMSD value of 0.29 nm. The drug irinotecan bound 2′OMTase demonstrated the average RMSD of 0.27 nm up to 40 ns, and thereafter the abrupt increase in RMSD was observed up to 50 ns. The 2′OMTas-teniposide complex displayed an average RMSD at 0.36 nm from 10 ns to 50 ns. The RMSD analysis confirmed that the apo and ligand-bound 2′OMTase attains a stable conformation throughout the MD simulation. Similarly, the ligand stability was also determined by estimating the RMSD throughout the simulation. The ligands viz. sinefungin, digitoxin, dihydroergotamine and teniposide showed RMSD value in the range of 0.2–0.3 nm confirming the stable conformation (Fig. 4B). However, the ligand, irinotecan, displayed the RMSD value in the range of 0.5–1.2 nm with a high degree of fluctuation, revealing that the ligand is flexible and unstable. The residue wise analysis of root means square fluctuation (RMSF) for apo and ligand-bound 2′OMTase based on Cα atoms revealed that the structural changes occur in the loop forming region and displayed the high flexibility (Fig. 4C). The protein regions displaying the high fluctuation are ranging from the amino acid 16–40, 72–80, 99–110, 132–147, 212–219, and 242–268 position shows the deviation up to 0.4 nm, apart from these regions other amino acids are showing variation from 0.1 to 0.25 nm. The ligand-bound complex viz. 2′OMTase-sinefungin, 2′OMTase-dihydroergotamine, 2′OMTase-irinotecan and 2′OMTase-teniposide, when compared with the apo-2′OMTase, displayed the high fluctuation in the temperature factor of the region mentioned above suggesting that these regions are mainly involved in the formation of catalytic cleft.

Fig. 4

Molecular dynamics simulation of Apo and ligand bound 2′OMTase A) Protein RMSD analysis, B) Ligand RMSD analysis, C) RMSF analysis, D) Radius of gyration analysis and E) Solvent Accessibility Surface Area Analysis.

The ligand-bound complex viz. 2’OMTase-sinefungin, 2’OMTase-digitoxin, 2’OMTase-dihydroergotamine, 2’OMTase-irinotecan and 2’OMTase-teniposide displayed the distance between the hydrogen bond forming atom was in the range of 0.28-0.33 nm. The estimated distance between the atoms follows the hydrogen bond formation rules and confirms the stable conformation of protein-ligand complex.

The radius of gyration analysis of apo 2′OMTase for the estimation of global compactness displayed the average Rg values of 1.90 nm up to 50 ns; however, a small increment in the Rg value was observed between 18 ns to 23 ns (Fig. 4D). The ligand-bound 2′OMTase complexes displayed the fluctuation in Rg value for an initial 3 ns, and thereafter the Rg value remains stable in the range of 1.85–1.89 nm throughout the simulation (Fig. 4D). The Rg analysis concludes that the protein remains in stable conformation resulted in achieving global compactness. Similarly, the estimated average solvent accessible surface area (SASA) for apo form of 2′OMTase up to 25 ns was 160 nm2 and thereafter, it reduces to 150 nm2 and remained unchanged throughout the MD simulation (Fig. 4E). The ligand-bound complex of 2′OMTase also displayed a similar type of pattern, and the average SASA was found to be in the range of 133–141 nm2.

The ligand-bound complex of 2′OMTase after the MD simulation was further subjected to the estimation of the interaction energy between protein and ligand molecule. The interaction energy of 2′OMTase for sinefungin was −54.35 kcal/mol, for digitoxin −69.47 kcal/mol, for dihydroergotamine −92.74 kcal/mol, irinotecan −75.34 kcal/mol, and teniposide −97.88 kcal/mol. The interaction energy calculation reveals that the digitoxin, dihydroergotamine, irinotecan, and teniposide has more affinity towards 2′OMTase and forms a stable conformation as compared to the reference ligand i.e. sinefungin.

Discussion

Methylation of mRNA Cap present at 5′ terminal is a crucial process, which protects its recognition as a foreign particle from the host immune system. 2′OMTase is a crucial enzyme involved in the methylation of 5′capped mRNA, which makes this enzyme a potential drug target for identification and screening of FDA approved drugs for repurposing. The modeled structure of the 2′OMTase protein displayed the conserved fold of the Class I MTase enzyme [12,13].

The quality assessment of the modeled structure by SAVES server reveals that the 100% non-glycine and non-proline residues attain the favorable dihedral, phi (ϕ), and psi (ψ) angles. The non-bonded interaction analysis of different types of atoms in the modeled structure when compared with the experimentally determined structure also suggested that the modeled structure is of excellent quality and can be used for further analysis. The multiple sequence alignment analysis of 2′OMTase with other homologous proteins revealed that the residues involved in catalytic tetrad formation i.e., Lys46, Asp130, Lys170, and Glu203, are strictly conserved in all the CoV strains [12,13]. Structure-based screening analysis of 3000 FDA approved drugs against 2′OMTase for repurposing helped in identifying the inhibitor molecules which can block the active site pocket of 2′OMTase enzyme. The initial top 22 FDA approved potential drug compounds included antiviral, alkaloid based drugs, cardiac glycoside, anticancer drugs, steroid-based drugs, and other drugs, etc. Five best compounds were selected based on the maximum number of conformations at a particular free energy of binding, ranging from −9.8 kcal/mol to −8.3 kcal/mol after re-docking analysis of the top 22 FDA approved drugs. These five best drug compounds were sinefungin, digitoxin, dihydroergotamine, irinotecan, and teniposide. These compounds interact through hydrogen bond and hydrophobic interactions predominantly with Asp130, Lys170, and Glu203, the amino acid residues involved in the formation of a catalytic tetrad of 2′OMTase. These amino acid residues viz. Asp130, Lys170, and Glu203 help in accommodating and stabilizing the ligands in the catalytic cleft. The apo and ligand-bound complexes were analyzed for conformational dynamics and stability by MD simulation studies. The RMSD analysis of apo and complex form revealed that except 2′OMTase-irinotecan complex, all the conformations attain the stable conformation. The RMSF analysis displayed that the regions of the protein ranging from 16–40, 72–80, 99–110, 132–147, 212–219, and 242–268 are mainly composed of β and γ turns and provide flexibility. The interaction energy analysis of the selected drugs revealed that digitoxin, dihydroergotamine, irinotecan, and teniposide have more affinity towards 2′OMTase and forms a stable conformation as compared to the earlier known inhibitor, i.e. sinefungin.

The nucleoside inhibitor, sinefungin, suppressed the replication of feline herpesvirus type 1 in the host feline kidney cells [17]. The alkaloid based drug, dihydroergotamine was not used as antiviral but for the migraine problems in HIV patients [18]. Cardiac glycosides such as digitoxin inhibit the human cytomegalovirus (HCMV) replication by inducing AMP-activated protein kinase (AMPK) activity and autophagy flux through the Na+, K+/ATPase α1 subunit activation [19]. Similarly, teniposide, the anticancer compound, inhibits the replication of SV40 by blocking the topoisomerase type II activity [20]. To date, all these potential FDA approved drugs are used for different applications and not used for the inhibition of 2′OMTase. The strong predicted binding of these approved drugs at the active site of 2′OMTase displayed that they can be used for inactivating the enzyme.

The SARS-CoV-2 infection has created a global pandemic and an urgent need for therapeutics. Repurposing of an existing drug based on the bioinformatics study could be a fast way to understand its antiviral efficacy against SARS-CoV-2. Although the understanding of SARS-CoV-2 molecular biology is limited, and there is a pressing need to design some potential therapeutics. The identification of some novel drugs that can display significant effect in silico can be further explored to validate its use clinically to control the disease and associated mortality. This study may pave the way in understanding the biology of SARS-CoV-2 and finding new possible drugs.

Declaration of competing interest

The authors do not have any potential conflict of interest.