Structural basis of SARS-CoV-2 3CLpro and anti-COVID-19 drug discovery from medicinal plants.

The recent pandemic of coronavirus disease 2019 (COVID-19) caused by SARS-CoV-2 has raised global health concerns. The viral 3-chymotrypsin-like cysteine protease (3CLpro) enzyme controls coronavirus replication and is essential for its life cycle. 3CLpro is a proven drug discovery target in the case of severe acute respiratory syndrome coronavirus (SARS-CoV) and Middle East respiratory syndrome coronavirus (MERS-CoV). Recent studies revealed that the genome sequence of SARS-CoV-2 is very similar to that of SARS-CoV. Therefore, herein, we analysed the 3CLpro sequence, constructed its 3D homology model, and screened it against a medicinal plant library containing 32,297 potential anti-viral phytochemicals/traditional Chinese medicinal compounds. Our analyses revealed that the top nine hits might serve as potential anti- SARS-CoV-2 lead molecules for further optimisation and drug development process to combat COVID-19.

Introduction

A novel coronavirus strain linked with fatal respiratory illness was reported in late 2019 [1]. Swift actions were taken by the Centre for Disease Control and Prevention (CDC), Chinese health authorities, and researchers. The World Health Organization (WHO) temporarily named this pathogen 2019 novel coronavirus (2019-nCoV) [2]. On January 10, 2020, the first whole-genome sequence of 2019-nCoV was released, which helped researchers to quickly identify the virus in patients using reverse-transcription polymerase chain reaction (RT-PCR) methods [3]. On January 21, the first article related to 2019-nCoV was published, which revealed that 2019-nCoV belongs to the beta-coronavirus group, sharing ancestry with bat coronavirus HKU9-1, similar to SARS-coronaviruses, and that despite sequence diversity its spike protein interacts strongly with the human ACE2 receptor [1]. On January 30, the WHO announced a Public Health Emergency of International Concern (PHEIC) for the 2019-nCoV outbreak. Later, the human-to-human transmission was confirmed. As of January 31, 51 whole-genome sequences of 2019-nCoV from different laboratories and regions had been submitted to GISAID database [4]. On February 12, the WHO permanently named the 2019-nCoV pathogen as SARS-CoV-2 and the causing disease as coronavirus disease 2019 (COVID-2019). The Chinese government’s swift actions helped them to control COVID-19 in China. However, SARS-CoV-2 affected several countries world-wide. On March 11, the WHO formally recognized the COVID-19 as a pandemic. By March 19, 2020, the global death toll reached 9,913, with 2,42,650 laboratory-confirmed cases. The case fatality rate among infected people varies in different countries. However, global case fatality rate is presently around 3.92% (calculated as deaths/[deaths + laboratory confirmed cases]).

Coronaviruses are single-stranded positive-sense RNA viruses that possess large viral RNA genomes [5]. Recent studies showed that SARS-CoV-2 has a similar genomic organization to that of other beta-coronaviruses, consisting of a 5′-untranslated region (UTR), a replicase complex (orf1ab) encoding non-structural proteins (nsps), a spike protein (S) gene, envelope protein (E) gene, a membrane protein (M) gene, a nucleocapsid protein (N) gene, 3′-UTR, and several unidentified non-structural open reading frames [3]. Although SARS-CoV-2 is classified into the beta-coronaviruses group, it is different from MERS-CoV and SARS-CoV. Recent studies highlighted that SARS-CoV-2 genes share <80% nucleotide identity and 89.10% nucleotide similarity with SARS-CoV genes [6,7]. Usually, beta-coronaviruses produce a ∼800 kDa polypeptide upon transcription of the genome. This polypeptide is proteolytically cleaved to generate various proteins. The proteolytic processing is mediated by papain-like protease (PLpro) and 3-chymotrypsin-like protease (3CLpro). The 3CLpro cleaves the polyprotein at 11 distinct sites to generate various non-structural proteins that are important for viral replication [8]. 3CLpro plays a critical role in the replication of virus particles and unlike structural/accessory protein-encoding genes, it is located at the 3′ end which exhibits excessive variability. Therefore, it is a potential target for anti-coronaviruses inhibitors screening [9]. Structure-based activity analyses and high-throughput studies have identified potential inhibitors for SARS-CoV and MERS-CoV 3CLpro [[10], [11], [12]]. Medicinal plants, especially those employed in traditional Chinese medicine, have attracted significant attention because they include bioactive compounds that could be used to develop formal drugs against several diseases with no or minimal side effects [13]. Therefore, the present study was conducted to gain structural insights into the SARS-CoV-2 3CLpro and to discover potent anti-COVID-19 natural compounds.

Materials and methods

Data collection

Whole-genome sequences of all SARS-CoV-2 isolates available till January 31, 2020, were downloaded from GISAID database (accession numbers and details are given in Table S1) [4]. The genome sequence of BetaCoV/Kanagawa/1/2020 (GISAID: EPI_ISL_402126) was incomplete, and the genome sequence of BetaCoV/bat/Yunnan/RaTG13/2013 (EPI_ISL_402131) was an old sequence (2013); therefore, these sequences were not included in our analyses. Gene sequences of 3CLpro were extracted from the whole-genome sequences and translated into protein sequences using the translate tool of the ExPASy server [14]. The first SARS-CoV-2 sequence (Wuhan-Hu-1; GSAID: EPI_ISL_402125) was used as a reference in our analysis.

Sequence analyses

In order to identify similar sequences and key/conserved residues, and to infer phylogeny, multiple sequence alignment of SARS-CoV-2 3CLpro followed by phylogenetic tree analyses were performed using T-Coffee [15] and the alignment figure was generated using ESPript3 [16]. Physicochemical parameters of SARS-CoV-2 3CLpro including isoelectric point, instability index, grand average of hydropathicity (GRAVY), and amino acid and atomic composition were investigated using the ProtParam tool of ExPASy [14].

Structural analyses

To probe the molecular architecture of SARS-CoV-2 3CLpro, comparative homology modelling was performed using Modeller v9.11 [17]. To select closely-related templates for modelling, PSI-BLAST was performed against all known structures in the protein databank (PDB) [18]. Chimera v1.8.1 [19] and PyMOL educational version [20] were used for initial quality estimation, energy minimisation, mutation analyses, and image processing.

Ligand database preparation and molecular docking

A comprehensive medicinal plant library containing 32,297 potential anti-viral phytochemicals and traditional Chinese medicinal compounds was generated from our previously collected data and studies [13,[21], [22], [23]], and screened against the predicted SARS-CoV-2 3CLpro structure. Molecular operating environment (MOE) [24] was used for molecular docking, ligand-protein interaction and drug likeness analyses. All analyses were performed using the same protocols that are already described in our previous studies [13,25,26]. The qualitative assessment of absorption, deposition, metabolism, excretion and toxicity (ADMET) profile of selected hits were predicted computationally by using ADMETsar server [27].

Molecular dynamics simulations

Explicit solvent molecular dynamics (MD) simulations were performed to verify docking results and to analyse the binding behaviour and stability of potential compounds using the predicted SARS-CoV-2 3CLpro homology model. GROMACS v5.1.4, GROMOS96 and the PRODRG server were employed to run 50 ns MD simulations [28,29] following the same protocol as described in our previous studies [13,30].

Results and discussion

Sequence and structural analyses

Multiple sequence alignment results revealed that 3CLpro was conserved, with 100% identity among all SARS-CoV-2 genomes. Next, the SARS-CoV-2 3CLpro protein sequence was compared with its closest homologs (Bat-CoV, SARS-CoV, MERS-CoV, Human-CoV and Bovine-CoV). The results revealed that SARS-CoV-2 3CLpro clusters with bat SARS-like coronaviruses and shares 99.02% sequence identity (Fig. 1A). Furthermore, it shares 96.08%, 87.00%, 90.00% and 90.00% sequence identity with SARS-CoV, MERS-CoV, Human-CoV and Bovine-CoV homologs, respectively (Fig. 1B). These findings were consistent with those of initial studies reporting that SARS-CoV-2 is more similar to SARS-CoV than to MERS-CoV, and shares a common ancestor with bat coronaviruses [1,3,31]. Analysis of physicochemical parameters revealed that the SARS-CoV-2 3CLpro polypeptide is 306 amino acids long with a molecular weight of 33,796.64 Da and a GRAVY score of −0.019, categorising the protein as a stable, hydrophilic molecule capable of establishing hydrogen bonds (Table 1).

Fig. 1

(A) Phylogenetic tree inferred from closest homologs of SARS-CoV-2 3CLpro. The maximum likelihood method was used to construct this tree. (B) Multiple sequence alignment of closest homologs of SARS-CoV-2 3CLpro sharing ≥70% sequence identity. (C) Cartoon representation of the SARS-CoV-2 3CLpro homodimer. Chain-A (protomer-A) is in multicolour and Chain-B (protomer-B) is in dark blue. The N-finger that plays an important role in dimerization maintaining the active conformation is shown in hot pink, domain I is coloured cyan, domain II is shown in green, and domain III is coloured yellow. The N- and C-termini are labelled. Residues of the catalytic dyad (Cys-145 and His-41) are highlighted in yellow and labelled. (D) Cartoon representation of the 3CLpro monomer model (chain/protomer-A) of SARS-CoV-2 superimposed with the SARS-CoV 3CLpro structure. The SARS-CoV 3CLpro template is coloured cyan, the SARS-CoV-2 3CLpro structure is coloured grey, and all identified mutations are highlighted in red. (E) Docking of 5,7,3′,4′-tetrahydroxy-2’-(3,3-dimethylallyl) isoflavone inside the receptor-binding site of SARS-CoV-2 3CLpro, showing hydrogen bonds with the catalytic dyad (Cys-145 and His-41). The 3CLpro structure is coloured dark blue, the 5,7,3′,4′-tetrahydroxy-2’-(3,3-dimethylallyl) isoflavone is orange, and hydrogen coloured maroon.

Table 1

Physicochemical parameters of SARS-CoV-2 3CLpro.

Parameters	SARS-CoV-2 3CLpro
Mol. Weight	33796.64 Dalton
No. of amino acids	306
Theoretical pI	5.95
Instability index (II)	27.65 (stable)
No. of negatively charged residues (Asp + Glu)	26
No. of positively charged residues (Arg + Lys)	22
Aliphatic index	82.12
Grand average of hydropathicity (GRAVY)	−0.019
Atomic composition	Carbon-1499; Hydrogen-2318; Nitrogen-402; Oxygen-445; Sulfur-22
Amino acid composition	Ala-17 (5.6%); Arg-11 (3.6%); Asn-21 (6.9%); Asp-17 (5.6%); Cys-12 (3.9%); Gln-14 (4.6%); Glu-9 (2.9%); Gly-26 (8.5%); His-7 (2.3%); Ile-11 (3.6%); Leu-29 (9.5%); Lys-11 (3.6%); Met-10 (3.3%); Phe-17 (5.6%); Pro-13 (4.2%); Ser-16 (5.2%); Thr-24 (7.8%); Trp-3 (1.0%); Tyr-11 (3.6%); Val-27 (8.8%); Pyl-0 (0.0%); Sec-0 (0.0%)

Next, for comparative modelling, BLAST [32] search identified SARS-CoV 3CLpro (PDB ID: 3M3V) as the best possible match in the PDB, with 100% query coverage, an E-value of 0.00, and 96.08% sequence identity. There were 12-point mutations (Val35Thr, Ser46Ala, Asn65Ser, Val86Leu, Lys88Arg, Ala94Ser, Phe134His, Asn180Lys, Val202Leu, Ser267Ala, Ser284Ala and Leu286Ala) between SARS-CoV and SARS-CoV-2 3CLpro enzymes (Fig. S1). Except for replacement of Leu with Ala at position 286, all other replacements conserve polarity and hydrophobicity. However, these mutations may affect 3CLpro structure and function. Therefore, the 3D structure of SARS-CoV-2 3CLpro was predicted. Firstly, a single chain monomeric model comprising all domains (Domain I = residues 8–100; Domain II = residues 101–183; Domain III = residues 200–303) was built (Fig. S2). N-terminal amino acids 1 to 7 form the N-finger that plays a significant role in dimerization and formation of the active site of 3CLpro. Domains I and II, collectively referred to as the N-terminal domain, include an antiparallel β-sheet structure with 13 β-strands. The binding site for the substrate is situated in a cleft between domains I and II. A loop from residues 184 to 199 joins the N-terminal domain and domain III, which is also called the C-terminal domain and comprises an anti-parallel cluster of five α-helices. The overall molecular architecture of SARS-CoV-2 3CLpro was in consistent with the crystal structure of SARS-CoV (PDB ID: 3M3V); the root mean square deviation (RMSD) between the homology model and the template was 0.629 Å. Structural and Ramachandran plot analyses revealed that 99% of residues are in favourable regions.

After quality assessment, individual chains were combined to form a homodimeric 3D structure, as shown in Fig. 1C. To facilitate other researchers, the predicted 3D model has been submitted to the Protein Model Database (PMDB) [33], and anyone can download/use the SARS-CoV-2 3CLpro final structure using PMDB ID: PM0082635. Furthermore, mutational analyses depicted none of the mutations affected the overall structure of SARS-CoV-2, which fully superimposed on the SARS-CoV 3CLpro structure (Fig. 1D). The results also revealed that SARS-CoV-2 has a Cys-His catalytic dyad (Cys-145 and His-41), consistent with that of SARS 3CLpro (Cys-145 and His-41), TGEV 3CLpro (Cys-144 and His-41) and HCoV 3CLpro (Cys-144 and His-41) [34]. These results revealed that the SARS-CoV-2 3CLpro receptor-binding pocket conformation resembles that of the SARS-CoV 3CLpro binding pocket and raises the possibility that inhibitors intended for SARS-CoV 3CLpro may also inhibit the activity of SARS-CoV-2 3CLpro.

Molecular docking

To test this hypothesis, we docked (R)-N-(4-(tert-butyl)phenyl)-N-(2-(tert-butylamino)-2-oxo-1-(pyridin-3-yl)ethyl)furan-2-carboxamide), a potential noncovalent inhibitor of SARS-CoV 3CLpro named ML188 [35], with the SARS-CoV-2 3CLpro homology model. We also docked ML188 with the SARS-CoV 3CLpro structure (PDB ID: 3M3V) as a reference, and ML188 bound strongly to the receptor binding site of SARS-CoV 3CLpro. The inhibitor targets the Cys-His catalytic dyad (Cys-145 and His-41) along with the other residues, and the docking score (S = −12.27) was relatively high. However, surprisingly, ML188 did not show significant binding to the catalytic dyad (Cys-145 and His-41) of SARS-CoV-2, and the docking score (S = −8.31) was considerably lower (Fig. S3). These results indicated that the 12-point mutations identified at previous step may disrupt important hydrogen bonds and alter the receptor binding site, thereby affecting its ability to bind with the SARS-CoV inhibitors.

Therefore, it is essential to discover novel compounds that may inhibit SARS-CoV-2 3CLpro and serve as potential anti-COVID-19 drug compounds. We developed a library from our previously published studies that contains numerous natural compounds possessing potential anti-viral activities and screened it against the SARS-CoV-2 3CLpro homology model. Recent drug repurposing studies proposed few drugs that target SARS-CoV-2 3CLpro and suggested that they could be used to treat COVID-19. Herein, we selected the best of these (Nelfinavir, Prulifloxacin and Colistin) from three different drug repurposing studies [36,37] and docked them as controls in the present study (Fig. S4). Our analyses identified nine novel non-toxic, druggable natural compounds that are predicted to bind with the receptor binding site and catalytic dyad (Cys-145 and His-41) of SARS-CoV-2 3CLpro (Table 2; Fig. S5). ADMET profiling of the selected hits is given in Table S2. Among these screened phytochemicals, 5,7,3′,4′-tetrahydroxy-2’-(3,3-dimethylallyl) isoflavone is an isoflavone extracted from Psorothamnus arborescens [38] that exhibited the highest binding affinity (−29.57 kcal/mol) and docking score (S = −16.35), and formed strong hydrogen bonds with the catalytic dyad residues (Cys-145 and His-41) as well as significant interactions with the receptor-binding residues Thr24, Thr25, Thr26, Cys44, Thr45, Ser46, Met49, Asn142, Gly143, His164, Glu166 and Gln189 (Fig. 1E). A literature review revealed that 5,7,3′,4′-tetrahydroxy-2’-(3,3-dimethylallyl) isoflavone has been successfully used as an anti-leishmanial agent [38], and it is also found in traditional Chinese medicine records [39]. Our screened phytochemicals displayed higher docking scores, stronger binding energies, and closer interactions with the conserved catalytic dyad residues (Cys-145 and His-41) than Nelfinavir, Prulifloxacin and Colistin. These results suggested that natural products identified in our study may prove more useful candidates for COVID-19 drug therapy.

Table 2

Summary of top ranked phytochemicals screened against SARS-CoV-2 3CLpro receptor binding site with their respective structures, docking score, binding affinity and interacting residues.

IDs	Phytochemical name	Plant source	Phytochemical structure	Docking score	Binding affinity (kcal/mol)	3CLpro residues interacting with phytochemical through H-bonding and other interactions
PubChem 11610052	5,7,3′,4′-Tetrahydroxy-2’-(3,3-dimethylallyl) isoflavone	Psorothamnus arborescens	Image 2	−16.35	−29.57	His41, Cys145, Thr24, Thr25, Thr26, Cys44, Thr45, Ser46, Met49, Asn142, Gly143, His164, Glu166, Gln189
PubChem 5281673	Myricitrin	Myrica cerifera	Image 3	−15.64	−22.13	Gly143, Cys145, His41, Thr24, Thr25, Thr26, Leu27, Cys44, Ser46, Met49, Leu141, Asn142, Ser144, His163, Glu166, Gln189
PubChem 6479915	Methyl rosmarinate	Hyptis atrorubens Poit	Image 4	−15.44	−20.62	Cys145, His41, Thr24, Thr25, Thr26, Cys44, Thr45, Met49, Leu141, Asn142, Gly143, Ser144, His163, His164, Glu166, Gln189
NPACT00105	3,5,7,3′,4′,5′-hexahydroxy flavanone-3-O-beta-D-glucopyranoside	Phaseolus vulgaris	Image 5	−14.42	−19.10	Met49, Cys145, His41, Thr24, Thr25, Thr26, Cys44, Ser46, Asn142, His164, Met165, Glu166, Gln189
PubChem 10930068	(2S)-Eriodictyol 7-O-(6″-O-galloyl)-beta-D-glucopyranoside	Phyllanthus emblica	Image 6	−14.41	−19.47	Thr24, Thr25, Gly143, Met49, Cys145, His41, Thr26, Cys44, Thr45, Glu166, Leu167, Gln189, Thr190, Ala191, Gln192
PubChem 5273567	Calceolarioside B	Fraxinus sieboldiana	Image 7	−14.36	−19.87	His41, Gly143, Cys145, Glu166, Thr24, Thr25, Thr26, Leu27, Ser46, Leu50, Leu141, Asn142, Ser144, His164, Met165, Gln189
PubChem 5318606	Myricetin 3-O-beta-D-glucopyranoside	Camellia sinensis	Image 8	−13.70	−18.42	Asn142, Glu166, Cys145, His41, Thr24, Thr25, Thr26, Thr45, Ser46, Met49, Leu141, Gly143, Ser144, His163, His164, Met165, Gln189
PubChem 11111496	Licoleafol	Glycyrrhiza uralensis	Image 9	−13.63	−19.64	Cys145, His41, Thr24, Thr25, Thr26, Cys44, Thr45, Met49, Gly143, His163, Met165, Glu166, Gln189
PubChem 6123095	Amaranthin	Amaranthus tricolor	Image 10	−12.67	−18.14	Thr26, Glu166, Cys145, His41, Thr24, Thr25, Cys44, Thr45, Asn142, His163, His164, Met165, Glu166, Leu167, Pro168, Gln189
PubChem 64143	Nelfinavir	Drugs used as control	Image 11	−12.20	−17.31	Met49, Met165, Glu166, Leu167, Pro168, Gly 170, Gln189, Thr190, Ala191
PubChem 65947	Prulifloxacin		Image 12	−11.32	−15.40	Gln192, Leu50, Met165, Glu166, Leu167, Pro168, Arg 188, Gln189, Thr190, Ala191
PubChem 5311054	Colistin		Image 13	−13.73	−18.57	Met49, Thr24, Thr25, Thr26, Thr45, Ser46, Glua47, Leu50, Asn142, Gly143, Met165, Glu166, Leu167, Pro168, Gln189, Thr190, Ala191, Gln192

∗3CLpro catalytic dyad (Cys-145 and His-41) residues are highlighted with bold font.

MD simulations

To further investigate the molecular docking results, the top three phytochemical complexes, namely 5,7,3′,4′-tetrahydroxy-2’-(3,3-dimethylallyl) isoflavone, myricitrin, and methyl rosmarinate, were subjected to 50 ns MD simulation. The root mean square deviation (RMSD), root mean square fluctuation (RMSF), radius of gyration (RoG) and hydrogen bond parameters were calculated. RMSD is an indicator of the stability of ligand-protein complexes. None of the complexes showed any obvious fluctuations, and all three were stable, with average RMSD values of 1.6 ± 0.02 Å, 1.5 ± 0.02 Å and 1.7 ± 0.02 Å for 5,7,3′,4′-tetrahydroxy-2’-(3,3-dimethylallyl) isoflavone, myricitrin, and methyl rosmarinate, respectively (Fig. 2A). RMSF is an indicator of residual flexibility. Minimal fluctuations were observed for myricitrin and methyl rosmarinate, and the overall complexes remained stable throughout the simulations. The functionally important catalytic dyad residues (Cys-145 and His-41) displayed stable behaviour, and fluctuations were observed toward the C-terminal end of the SARS-CoV-2 3CLpro molecule (Fig. 2B). RoG is an indicator of protein compactness, stability, and folding, and the results suggested normal behaviour for all three complexes; all remained compact and stable throughout the 50 ns simulations (Fig. 2C). In addition, hydrogen bonds, which are the main stabilising interactions factors in proteins, suggested that the SARS-CoV-2 3CLpro internal hydrogen bonds remain stable throughout the simulation, with no obvious fluctuations (Fig. 2D). These results confirmed our findings and further indicated that these compounds may serve as potential anti-COVID-19 drug sources.

Fig. 2

(A) Root mean square deviation (RMSD), (B) root mean square fluctuation (RMSF), (C) potential energy and (D) Hydrogen Bond interactions for all three complexes over the 50 ns simulation.

Conclusion

In conclusion, our study revealed that 3CLpro is conserved in SARS-CoV-2. It is highly similar to bat SARS-like coronavirus 3CLpro, with some differences from other beta-coronaviruses. We predicted the 3D structure of the SARS-CoV-2 3CLpro enzyme, and the findings may help researchers working on COVID-19 drug discovery. Despite significant overall similarity with the SARS-CoV 3CLpro structure, the SARS-CoV-2 3CLpro substrate binding site had some key differences, which highlighted the need for rapid drug discovery to address the alarming COVID-19 pandemic. Medicinal plant compounds have already been used to successfully treat numerous viral diseases. Herein, we screened a medicinal plant database containing 32,297 potential anti-viral phytochemicals and selected the top nine hits that may inhibit SARS-CoV-2 3CLpro activity and hence virus replication. Further in-vitro and in-vivo analyses are required to transform these potential inhibitors into clinical drugs. We anticipate that the insights gained in the present study may prove valuable for exploring and developing novel natural anti-COVID-19 therapeutic agents in the future.

Conflicts of interest

The authors declare that there are no conflicts of interest.