Computational investigation of phytoalexins as potential antiviral RAP-1 and RAP-2 (Replication Associated Proteins) inhibitor for the management of cucumber mosaic virus (CMV): a molecular modeling, in silico docking and MM-GBSA study

Abstract The Replication Associated Proteins (RAP-1 and RAP-2) encoded by CMV ORF 1a and ORF 2a are required for the different stages of the viral replication cycle; being multi-functional, they are good inhibitory targets for anti-CMV compounds. As a new perspective for sustainable crop improvement, we investigated the natural plant-based antimicrobial phytoalexins for their anti-CMV potential. Here, we modeled and predicted the functional domains of RAP-1 and RAP-2, docked with a ligand library comprising 128 phytoalexins reported with broad-spectrum activity, determined their binding energies (BEs), molecular interactions, and inhibition constant (Ki), and compared with the reference plant antiviral compounds ribavirin, ningnanmycin, and benzothiadiazole (BTH). Further, the change in Gibb’s free energy of binding (ΔG) and the per residue contribution of the selected top-scored ligand molecules was assessed by the prime MM-GBSA approach. Our results revealed RAP-1 as a discontinuous two-domain and RAP-2 as a multi-domain protein. The compounds glyceollidin (9.8 kcal/mol) and moracin D (7.8 kcal/mol) topped the list for RAP-1 and RAP-2 protein targets respectively and also, the lead molecules had energetically more favorable and comparative ΔG values than the top-scored plant antiviral agent ningnanmycin. The evaluation of in vitro toxicity and agrochemical-like properties showed the least toxicity of these anti-CMV compounds. Taken together, our results provide new insights in understanding the inhibitory effects of phytoalexins towards the RAP proteins and could be employed as new promising anti-CMV candidate compounds for their application in agriculture as biopesticides to combat the CMV disease incidence. Communicated by Ramaswamy H. Sarma The CMV replication cycle and the proposed effect of viral replication inhibitors. Upon entry, CMV un-coats the coat protein (CP), co-localizes in the endoplasmic reticulum for its replication, mRNA synthesis, and translation. The viral RNA particles are assembled with the CP in the cytoplasm and released. The proposed anti-CMV strategy of phytoalexins as potential inhibitors for targeting the two viral enzymes RAP-1 and RAP-2 may likely inhibit the key steps 3, 4, 5, and 6 essential for the successful completion of the CMV life cycle and minimize the viral spread.


Introduction
Agricultural crop production is threatened by several abiotic and biotic stresses and their combinations thereof. It is estimated that among the biotic stresses, about 40% of total crop losses are accounted for by viral infection. The majority of viruses that infect agricultural plants (at least 450 different species) are RNA viruses . Cucumber Mosaic Virus (CMV) (genus: Cucumovirus, family: Bromoviridae and alphavirus-like superfamily) is identified as one of the most devastating plant RNA viruses with an unusually broad host range covering 1200 crop species and 100 families. The crop losses accounted for the CMV infection are enormously making it a highly attractive virus with a significant impact on agriculture crop production. The genome is a tripartite positive single-stranded RNA where RNA 1 and RNA 2 ORFs encode Replication Associated Protein 1 (RAP-1) and Replication Associated Protein 2 (RAP-2) subunits; their compatible interactions along with the host factors are very essential for the in vivo and in vitro viral replication cycle, transcription, and translation (Tzean et al., 2019). However, the non-structural RAP proteins are largely uncharacterized and are rarely targeted for CMV disease management (Sal anki et al., 2018).
Although considerable progress is being made for crop improvement, there are no direct methods to combat viruses, unlike bacterial and fungal pathogens. Many conventional methods are employed to combat the CMV infection; cultural practices, vector eradication using synthetic pesticides which are not long term sustainable, development of CMV resistant cultivars but are limited due to lack of complete resistance to CMV in many crops. The quanta of compounds identified as anti-CMV are seldom rare due to the lack of molecular mechanisms underlying them. Some of the commercially available plant antiviral agents include ribavirin, ningnanmycin, and benzothiadiazole (BTH), but are not large-scale field-deployable due to the cost and safety concerns (Guo et al., 2021;Song et al., 2014). The modern approaches to plant protection via genetic engineering and CRISPR genome editing also demand identifying several plant-based candidate compounds with promising antiviral activity.
In this context, we selected phytoalexins, which are versatile, low molecular weight induced plant-based antimicrobial compounds reported with a wide spectrum of biological activities (Jeandet, 2015). They are an integral component of the biochemical defense arsenal of the plants; known to exhibit biocidal as well as biostatic action on several pathogens including fungal, bacterial, viruses, nematodes, arachnids, other animals, and plants, and their role in disease resistance is also well documented (Guest, 2017). Several research findings have indicated that phytoalexins have plant antiviral potential, for example; accumulation of phytoalexin glutinosone in Nicotiana glutinosa upon TMV infection (Burden et al., 1975), capsidiol accumulation in Nicotinia clvelandi and N tabacum during Tobacco necrosis virus infection (TNV) (Chaube and Pundhir 2005), trans-resveratrol accumulation in grapevine leafroll-associated virus 3 (GLRaV-3)-infected plant leaves (Bona et al., 2020), but the underlying mechanism has remained unknown. There have also been few attempts to determine the suppression of viral replication by phytoalexins. Studies conducted by Hammerschmidt (1999) and Li et al. (2014) have revealed the ability of glyceollin to restrict TNV multiplication in soybean leaves and gossypol Schiff base functional analogs exhibiting high inhibitory activities towards TMV viral replication.
Computational approaches can aid in understanding the molecular interactions of anti-CMV compounds with the viral target proteins. In this study, we investigated the molecular interactions of the selected phytoalexins towards RAP-1 and RAP-2 protein targets through a molecular docking and MM-GBSA approach to determine their potential inhibitory effects. Our goal of the current research was to identify the anti-CMV potential of phytoalexins for their use as an alternative, natural and safe biopesticides to combat viral diseases affecting crop plants.

Physicochemical parameters, secondary structure, and disulfide bond (S-S) prediction
The RAP-1 comprising 993 (NCBI Accession number: NP_049323.1) and RAP-2 protein with 857 amino acids (NCBI Accession Number: NP_049324.1) were retrieved from the NCBI database (https://www.ncbi.nlm.nih.gov/). The primary sequence analyses of both the proteins were calculated using protein identification and analysis tool ExPASy Protparam server (https://web.expasy.org/protparam/). The computed physicochemical parameters include molecular weight (Mw), theoretical isoelectric point (pI), amino acid composition, extinction coefficient, instability index (II), aliphatic index (AI), and grand average of hydropathicity (GRAVY) (Gasteiger et al., 2005) (https://web.expasy.org/protparam/). The secondary structure was predicted using the Self-Optimized Prediction Method with Alignment (SOPMA) (https://npsa-prabi.ibcp.fr/cgi-bin/npsa_automat.pl?page=/ NPSA/npsa_sopma.html). The SOPMA predicts the secondary structure by consensus prediction from multiple alignments of the set of proteins belonging to the same family (Geourjon & Deleage, 1995). The disulfide bonds for both the target proteins were predicted using DiANNA 1.1 webserver (Ferr e and Clote 2006) which determines the cysteine oxidation state and disulfide connectivity topology of the input amino acid sequence based on the implementation of the neural network of Fariselli et al. (1999).

Quality assessment and model validation
3D structures of RAP-1 and RAP-2 proteins were predicted by Iterative Threading Assembly Refinement (I-TASSER, version 5.1) (https://zhanglab.ccmb.med.umich.edu/I-TASSER/). It is an automated server and a hierarchical protein structure modeling approach based on the statistical significance of the profile-profile alignment (PPA) threading alignments and the structure convergence of the Monte Carlo simulations Zhang & Skolnick, 2004). In brief, the output file comprises a predicted confidence score (C-score) for each of the five models in the range of À2 to þ5, which predicts the reliability and the quality of the model, with the RMSD and TM score for the first model Zhang, 2008). The models with the lowest confidence score (C-score) were subsequently selected for the structural refinement using an Ab initio loop and terminus-based GalaxyRefine 2 webserver (http://galaxy.seoklab.org/cgi-bin/ submit.cgi?type=REFINE2) (Ko et al., 2012). The quality of the final model was assessed by the PROCHECK program of Structural Analysis and Verification (SAVES 6.0) (https://saves. mbi.ucla.edu/) tool to determine the stereochemical quality of the model by generating a Ramachandran plot as previously described (Laskowski et al., 1993;Sidhu et al., 2020;Sivaramakrishnan et al., 2012).

Preparation of receptors and ligands
Both the receptors and ligands were prepared in the PDBQT format using ADT tools (ADT) (v.1.5.6) of the MGL software package. A semi-flexible docking protocol was adopted where the receptors were kept rigid and the ligands were kept flexible (Gowthaman et al., 2008). Briefly, protein receptors were prepared by removing water molecules and hydrogen bonds, merging non-polar atoms, adding Kollman and Gasteiger charges. On the other hand, the selections of ligands were based on an extensive literature survey and a ligand library comprising 128 phytoalexins was prepared. The compounds ribavirin, BTH, and ningnanmycin reported as commercial plant antiviral agents in the scientific literature were chosen as reference compounds (Supporting information Table S1). All the ligand molecules were optimized with the following procedure; the 3D structures were retrieved from the Pubchem ligand database (https://pubchem.ncbi. nlm.nih.gov/) and saved in Structured Data Format (SDF). The conversions of SDF to Protein data bank (PDB) format files were done using Open Babel software tools (O'Boyle et al. 2011); detected the number of rotatable bonds plus torsional degrees of freedom (TORSDOF) and converted into PDBQT file format using ADT tools. A total of 128 configuration files were generated along with the reference ligand compounds to perform an independent docking procedure with AutoDock Tools (ADT) (v.1.5.6) (Sanner 1999).

Molecular Docking
AutoDock tools were used for preparing, running, and analyzing the docking of the target proteins RAP-1 and RAP-2 with the prepared ligand library (http://vina.scrip ps.edu/) (Trott & Olson, 2010). AutoDock requires pre-calculated grid maps, one for each atom type; present in the ligand being docked and it stores the potential energy arising from the interaction with the macromolecule. This grid must surround the region of interest in the macromolecule. Both the models had a discrete cubic grid box with dimensions of 28 Â 28 Â 28 points (XYZ dimensions) and a resolution of 1 Å spacing generated using the 'grid box generation platform' of the MGL tools, large enough to allow for the free rotation of the ligand molecules around the selected binding site residues (Natesh et al., 2021). The X center, Y center, and Z center were 101.74, 70.901, and 87.263 for RAP-1 and 98.96, 82.096, and 113.323 for RAP-2 protein respectively. During the docking process, a maximum of 10 conformers was considered for each compound. The docking protocol was repeated at least twice to ensure the accuracy of the generated results.

Ligand-receptor interaction and inhibition constant calculation
Protein-ligand complexes were visualized and analyzed using both PyMol (The PyMOL Molecular Graphics System, version 2.4.1 Schrodinger, LLC) and Discovery studio visualizer (v21.1.0.20298). The protein-ligand complexes were evaluated based on the binding affinity and interaction of the ligand with the receptor. Further, the inhibition constant (Ki) was calculated according to the method of Onawole et al. (2018) and the equation is as follows, Ki ¼ 10 ðBinding energy=1:366Þ

Molecular Mechanics -Generalized born model and solvent accessibility (MM-GBSA) energy calculations
The molecular docking at Standard precision (SP) was carried out using the Glide module, and MM-GBSA calculations were carried out using the Prime module of Schrodinger suites.
The hardware used for this study was from DELL brand with Intel Core i7-6700 clocked at 3.40 GHz X 8 processor with 8 GB RAM running on 64-bit CENT OS Linux operating system.

Protein preparation
The three-dimensional (3D) structure of the proteins RAP-1 and RAP-2 were imported in the Protein preparation wizard in Maestro. The protein was pre-processed by performing assign bond orders, removing original hydrogens, adding hydrogens, creating zero-order bonds to metals, creating disulfide bonds, converting seleno-methionines to methionines, deleting water molecules beyond 5 Å from hetero groups, and generating het states using Epik at pH: 7.0 p/-2.0. Then, the protein was refined by performing H-bond assignment by sample water orientations, using PROPKA pH: 7.0, and structure optimization was carried out. The restrained protein minimization was performed using the OPLS3 force field (Harder et al., 2016;Rawat & Verma 2020a, b).

Receptor grid generation
A receptor grid of size 25 Å was generated after preparing the proteins at the centroid of the active site.

Preparation of ligand
All the molecules were built, optimized as a 3D structure, and energy minimized using ligprep and confgen function in Maestro. The library was saved in .maegz format.

Molecular docking
Subsequently, molecular docking was performed using the Glide at Standard Precision (SP) algorithm to evaluate the protein-ligand interaction and their binding affinities (Friesner et al., 2004;Kant et al. 2020). For this study, the add Epik state penalties to docking score and enhanced planarity of conjugated pi group functions were selected.

MM-GBSA calculations
The free energies of binding for all complexes were calculated using the Molecular Mechanics: Generalized Born model and Solvent Accessibility (MM-GBSA) calculations using the Prime module. The change in binding energy (DG) Bind calculation was carried out for top hits and standards and per residue interaction, energy plots for these ligands were plotted, the changes in the free energies of the protein-ligand complexes (DG Bind) were calculated using the MM-GBSA Prime module of the Schr€ odinger suite 2018-4. The calculated DG Bind, per residue contribution in the total binding interaction energies and the energy plots for the top-ranked ligands of RAP-1 and RAP-2 protein targets, were compared with the reference compound ningnanmycin. The Prime MM-GBSA uses VSGB2.0 implicit solvent model and the OPLS2005 force field to compute the bound and unbound molecules involved in the binding process (Sirin et al., 2014;Rawat & Verma 2020a, b).

Toxicity prediction
The application of phytoalexins in green agriculture as anti-CMV compounds requires regular evaluation and a thorough food and feed safety assessment to determine their impacts on humans and environmental health. Computational approaches using in silico studies offer a simple and better alternative solution to the time-consuming experimental approaches. The in vitro toxicity prediction for the top-ranked phytoalexins including the median lethal dose (LD 50 ), toxicological endpoints, and targeting pathways was determined using the ProTox II server (https://tox-new.charite.de/protox_ II/) (Banerjee et al., 2018).

Prediction of physicochemical properties of ligands
Physicochemical properties of top-ranked phytoalexins were calculated using MarvinSketch (v 21.4). Physicochemical descriptors such as logP value, donor sites, acceptor sites, 2D Polar Surface Area (PSA), molar refractivity, Van der Waals 3D surface area, and water solubility were predicted by

Results and discussion
3.1. Physicochemical, secondary structure, and S-S bond characterization The primary sequence analyses computed by Protparam and the comparisons of the estimated parameters for both the viral replicase RAP-1 and RAP-2 proteins are summarized in Table 1. The Mw of RAP-1 (111.43 kDa) was larger than the RAP-2 protein (96.71 kDa). Comparative analyses showed significant differences in their theoretical pI and amino acid compositions; RAP-1 had a basic pI of 7.28 and RAP-2 with acidic pI 5.32, respectively. The variations in the pI values must reflect the changes in the amino acid composition. The aliphatic non-polar amino acid residues Alanine (Ala), Valine (Val), Leucine (Leu); and polar uncharged residues Serine (Ser) and Threonine (Thr) were in higher proportions for RAP-1 protein whereas amino acid residues Serine (Ser), Leucine (Leu) and negatively charged, acidic amino acid residue Aspartic acid (Asp) was in abundance for RAP-2 protein (Supporting information Figure S1). The extinction coefficient value was significantly higher in RAP-1 than in RAP-2 protein.

Domain boundary locations
The domain patterns and domain boundary locations of RAP-1 and RAP-2 proteins were determined by FUpred algorithm. Briefly, it generates a Folding Unit score (FUscore) which maximizes the number of intra-domain contacts and minimizes the number of inter-domain contacts. The domain boundaries are predicted in the regions with low FU scores for both continuous and discontinuous proteins (Zheng et al., 2020). Based on the prediction results, RAP-1 was a discontinuous two-domain protein with domain boundary locations at (D1-1)1-601, (D1-2) 637-993, and (D2) 602-636 respectively (Figure 1 a-c). These results were comparable with the previous findings suggesting two conserved domains in the members of the alphavirus-like superfamily; N-terminal with putative RNA capping-like domain demonstrated with m7G methyltransferase and covalent GTP binding (Ahola & Ahlquist, 1999;Kong et al., 1999) and a Cterminal with putative RNA helicase-like domain (Gorbalenya et al., 1988). Likewise, RAP-2 was predicted as a multi-domain discontinuous protein with 7 domain boundary locations at  (Argos, 1988;Haseloff et al., 1984).

3D Structural prediction and refinement
The 3D structures were determined by the I-TASSER server. It generates structural conformations, then uses the SPICKER clustering program to cluster all structures based on the pairwise structural similarity (Zhang, 2008). The model with the lowest C-score value relative to the native was selected as the best model. The C-scores, TM scores, and RMSD values of RAP-1 and RAP-2 proteins were; À0.06 and À0.89, 0.71 ± 0.12 and 0.06 ± 0.14, 9.0 ± 4.6 Å, and 10.7 ± 4.6 Å respectively. The top predicted models of RAP-1 and RAP-2 proteins were submitted separately to the GalaxyRefine 2 server for refinement. It rebuilds the side-chain and performs side-chain repacking and structure relaxation by molecular dynamic simulation (Ko et al., 2012). After the refinement process, refined models were submitted to SAVES 6.0 tool for validation. Most of the amino acid residues of RAP-1 and RAP-2 proteins were in the favored and allowed regions (95.8% and 96.1%) with 2.2% and 2.5% residues in the disallowed region respectively indicating a good quality of the refined models (Figure 2a, b). Moreover, both the proteins had the majority of the amino acid residues favoring b sheets and right-hand a helix. The cartoon and surface representation of RAP-1 and RAP-2 proteins along with their predicted binding sites are illustrated in Figure 3.

Molecular interaction of phytoalexins with RAP-1 and RAP-2 protein targets
The RAP-1 and RAP-2 proteins of CMV are known to have interconnected functions and they are thus important candidate targets and discovering inhibitor compounds are essential to disrupt and block the key steps involved in the viral replication life cycle. In the present molecular docking study, we docked RAP-1 and RAP-2 with 128 phytoalexins belonging to 28 families shown to display broad-spectrum activity to identify prospective anti-CMV inhibitors. In addition to the standard plant antiviral drugs, a total of 128 natural plant-based phytoalexins were docked with both the protein complexes. The BE scores of the top-ranked phytoalexins for the RAP target proteins with their potential interacting residues are summarized in Supporting information Table S3 and S4 respectively. The BEs of the ligands ranged from À5.4 to À9.9 kcal/mol for RAP-1 and 4.7 to À7.8 kcal/ mol for RAP-2 protein (Supporting information Table S5). Among the standard plant antiviral compounds, ningnanmycin, a cytosine nucleoside type antibiotic, isolated from Streptomyces noursei var. x ichangensis, showed the highest binding affinity (BE: À7.1 kcal/mol for RAP-1 and À6.2 kcal/ mol for RAP-2 target) than other standards. Moreover, from our results, nearly 87 and 72 phytoalexins had encouraging binding affinities for the RAP targets than the top-scored antiviral compound ningnanmycin (Supporting information  Table S5). Higher the binding affinities towards the target protein better the possibility of it as a more potent antiviral compound. The comparisons of the interacting residues of the top-ranked phytoalexins for the RAP targets along with the standard antiviral compound ningnanmycin are illustrated in Figures 4 and 5. Glyceollidin; a pterocarpan isoflavonoid phytoalexin originally isolated from Glycine max (Lygin et al., 2010) had the strongest binding affinity (BE: À9.9 kcal/mol) characterized to form two H-bond with Thr520, Lys523, and several hydrophobic interactions with amino acid residues; Asp521, Gly522, Ile524, Leu525 Ala526, Ala527, Thr554, Pro555, Thr558, Glu784, Gly809, Asp810, Gln813, Gly937, Ile938 (Supporting information Table S3, Figure 4b). Similarly, pisatin, an isoflavonoid phytoalexin originally isolated from Pisum sativum (Parniske et al., 1991) had higher binding energy scores (BE of À9.7 kcal/ mol) forming a single H-bond with Ser939 residue and nine non-covalent interactions with Lys523, Ile524, Leu525, Ala526, Ala527, Thr558, Glu784, Arg843, Lys970 amino acid residues (Supporting information Table S3, Figure 4c). The compounds glyceollidin, pisatin, maackiain, and phaseollin had Lys523 as a critical amino acid residue forming either H-bonded or hydrophobic interactions. In addition to the enzymatic function, RAP-1 also influences viral systemic movement and symptom severity (Watt et al., 2020). Though only limited attempts have been made for the plant antiviral potential of phytoalexins, their accumulation in response to viral infection has been reported by earlier studies. For example, phaseollin and pisatin accumulation in Phaseolus vulgaris and Pisum sativum in response to tobacco necrosis virus (TNV) supports our research findings (Bailey, 1973). Therefore, binding of these compounds to the RAP-1 could alter its function and may interrupt it from binding with the RAP-2 and likely block the virus multiplication and restrict its movement and minimize symptom severity. From our docking results, it was evident that all the top-ranked compounds for RAP-1 were isoflavonoid phytoalexins of the Fabaceae family and several research studies have demonstrated their significance in plant disease resistance (Jeandet, 2015).  Similarly, among all the docked compounds, Moracin D, a furanopterocarpan isolated from the Moraceae family (Takasugi et al., 1979) was characterized with the lowest binding energy showing the strongest binding interaction with RAP-2 protein (BE: À7.8 kcal/mol) forming three conventional H-bonds with Phe476, Asp610 amino acid residues, and several hydrophobic bonds with Gln471, Pro474, Leu475, Phe476, Thr477, Tyr558, Ser608, Gly609, Asp611, Lys648, Leu650, Lys674 amino acid residues (Supporting information Table S4, Figure 5b). Likewise, sanguinarine, an alkaloid phytoalexin isolated from Papaver somniferum had the second-highest BE score value (-7.4 kcal/mol) forming two H- bonds with Arg560, Lys648 and several hydrophobic interactions with Ile453, Gln471, Pro474, Leu475, Phe476, Tyr558, Arg560, Asp561, Ser608, Gly609, Asp610, Lys648, Leu650, Lys674 amino acids respectively (Supporting information Table S4, Figure 5d). The RAP-2 contains the domains specifying putative RNA-dependent RNA polymerase (RdRp) activity. RdRp of the RNA viruses is required for all the steps of viral replication. These enzymes lack proofreading activity and are highly error-prone (Smith, 2017), with high mutation rates and in such instances, the use of single antiviral compounds such as antibiotics will be ineffective. A combination of phytoalexins with broad-spectrum activity would be an ideal choice for the inhibition of RdRp to avoid the resistance development due to mutation. Among all the known RdRps, the greatest structural conservation is shown in the palm domain, comprising one of the five motifs; the Mg21-binding GDD motif, which is essential for the catalytic activity (Gorbalenya et al., 2002). Here, in our case, it was identified in domain 6 with positions at 609-611 proceeded by a glycine residue, and all the top-ranked phytoalexins formed either covalent or non-covalent interaction with this motif suggesting their potential ability to block the catalytic activity of RdRp and thus leading to the loss of its function. The chemical structures of the top-scored ligands for the RAP targets are given in Figure 6. It is also noteworthy to mention that, sixteen phytoalexins exhibited higher binding affinities towards both the RAP target proteins (Supporting information Table S6). The comparison of binding energies (BEs) of reference compound ningnanmycin with top-ranked phytoalexins for RAP proteins along with their corresponding counterpart BE values are given in Supporting information Figure  S2. Strong binding of these compounds either towards RAP-1 or RAP-2 or both proteins could influence their compatible interactions and may likely halt the CMV replication and block the viral spread.
The Ki values calculated for the top-ranked phytoalexins of both the targets along with the standard plant antiviral compound ningnanmycin are given in Supporting    information Tables S3 and S4. Ki is the concentration required to produce half the maximum inhibition and 1-40 mM is considered as a preferred range for a hit compound (Naidoo et al., 2020). The lower the Ki value, the higher the potency of the compound to inhibit the target molecule. Glyceollidin had the lowest Ki of 0.0563 lM for RAP-1 and Moracin-D with 1.94 lM for RAP-2 target. It is worth noting that, all the top-ranked compounds had Ki values lower than ningnanmycin, either in sub-nanomolar or micromolar concentrations, qualifying them as hit compounds with higher potency to inhibit the target molecules.
3.5. Calculation of change in Gibb's free energy of binding and per residue energy contribution using molecular Mechanics -Generalized born and surface area continuum solvation (MM-GBSA) To compare the binding of various molecules, the DG Bind for each ligand-receptor complex was calculated using the Molecular Mechanics-Generalized Born model and Solvent Accessibility (MM-GBSA) energy calculations which utilizes implicit solvation and better predicts the binding affinity of the ligands than molecular docking studies. The DG Bind for  each ligand with proteins RAP-1 and RAP-2 are mentioned in Tables 2 and 3, respectively. To our delight, the data revealed 2 hits namely glyceollidin and pisatin displaying better DG Bind values for protein RAP-1 than the reference antiviral compound, ningnanmycin. The DG Bind against RAP-1 for the molecules, glyceollidin, pisatin, and ningnanmycin were À43.02, À41.46, and À35.90 kcal/mol, respectively (Table 2 and Figure 7). On the other hand, for RAP-2 the molecule piceid attained a comparable DG Bind as the reference, ningnanmycin. The DG Bind against RAP-2 for the molecules, piceid and ningnanmycin were À41.82 and À41.13 kcal/mol, respectively (Table 3 and Figure 8). The detailed 3D interaction plots of the top-scored ligand molecules with the RAP-1 and RAP-2 targets are given in Supporting information Figures S3 and S4. Further, to understand the energy contributions of individual residues inside the active site, the per residue interaction energies were calculated for all the complexes. The amino acids which were involved in interactions with the inhibitors were contributing much towards the binding energy. In the case of RAP-1 complexes, the major contributing residues include Thr520, Gly522, Lys523, Thr558, Asp783, Glu784, Gly809, and Gln813. Among all particularly the energy contributions from Lys523 were considerably higher reflecting a strong H-bond interaction in the case of glyceollidin and ningnanmycin; and pi-stacking interaction in the case of phaseollin. Also, the residue Glu784 own significant binding contribution as an H-bond interacting residue; and Thr520, Gly522, Asp783, Gly809, and Ser939 were other residues inside the cavity that contributed through H-bonding interactions (Table 2 and Figure 9).
In the case of the RAP-2 complexes, the active site was filled with residues such as Gln533, Ile536, Phe556, Tyr558, Arg560, Asp561, Ile645, and Lys648. From the interaction diagrams, it was observed that most of these residues participated with H-bonding interactions. The molecule of interest, piceid acquired its strength from the contribution of double H-bonding interaction with residue Asp561, one Hbonding interaction with Lys648, and one pi-stacking interaction with residue Phe556 (Table 3). From the energy profile of contributing residues, it was observed that the highest share is from Asp561 and Lys648 (Figure 10).

Toxicity prediction
The in vitro toxicity levels of the phytoalexins were predicted with the ProTox-II server. ProTox-II incorporates molecular similarity, pharmacophores, fragment propensities, and machine-learning models for the prediction of various toxicity endpoints; such as acute toxicity, hepatotoxicity, cytotoxicity, carcinogenicity, mutagenicity, immunotoxicity, adverse outcomes pathways (Tox21), and toxicity targets (Banerjee et al., 2018). The in vitro toxicity prediction of the top-ranked phytoalexins showed LD 50 (mg/kg) ranging from 500 À 5010, falling under three different classes of toxicity (IV, V, and VI respectively). The top-ranked compound for RAP-1 protein, glyceollidin had LD 50 value 2500 (mg/kg) falling under class V toxicity with no hepatotoxicity, mutagenicity, carcinogenicity, and cytotoxicity except being immunotoxic. Similarly, Moracin D a top-ranked compound for RAP-2 protein had LD 50 of 5000 (mg/kg) falling under class V toxicity and immunotoxic. All the compounds had no hepatotoxicity and cytotoxicity and the majority were nonmutagenic indicating their least toxicity for human consumption and environment, than the commercially available synthetic pesticides (Table 4). The adverse outcomes pathways (Tox21) and toxicity targets for the top-ranked phytoalexins are listed in Supporting information Table S7.

Evaluation of physicochemical properties and analysis of agrochemical-like activity
Chemaxon tool was employed to evaluate the physicochemical properties of the top identified molecules exhibiting antiviral properties. The descriptors such as lipophilicity (expressed as consensus log P), molecular weight (Mw), donor sites, acceptor sites, 2D polar surface area (PSA), molar refractivity, Van der Waals 3D surface area, and water solubility were calculated to know the chemical nature of the bioactive compounds. Lipinski's rule of five (RO5) is considered as the reference for defining physicochemical and structural   parameters for oral drug bioavailability which is also applicable for evaluating the agrochemical-like property of any bioactive agent (Sidhu et al., 2020). As per rule, the logP value should be less than 5, the H-bond donor should be less than 5, the H-bond acceptor should be less than 10 and molecular weight should be less than 500 (Lipinski et al., 1997). The main aim of determining the physicochemical properties of the phytoalexins is to evaluate their potential to be applied as biopesticides in agriculture. The molecular weight of the top identified compounds ranged from 268.22-390.4 respectively. All the compounds exhibited a polar surface area value less than 140 Å 2 indicating good cell membrane permeability (Daina et al., 2017). Moreover, the molar refractivity values for the top-scoring phytoalexins were in the range of 69.2-99.6 displaying better agrochemical molecules-like properties (Pathak et al., 2017). Solubility of a molecule indicates the ease of handling and formulation of the phytoalexins as a biopesticide for agrochemical purposes. The compound piceid had the highest solubility in water whereas, other compounds were either moderately or less soluble in water. The list of physicochemical parameters calculated for all the top identified phytoalexins is given in Table 5.

Conclusion
The excessive use of synthetic agrochemicals in agriculture to improve crop productivity has resulted in serious consequences such as climate change which is one of the global environmental concerns, and also has led to other effects namely soil degradation, pesticides residual effect, insect resistance, and agro-ecological imbalance. It is therefore an alarming situation and necessitates prospecting natural plant-based products with antimicrobial activity to combat various stresses. Plants have inbuilt defense mechanisms naturally to overcome the invading pathogens, in this regard we explored the possibility of using plant-based induced biochemical defense compounds; phytoalexins reported with broad-spectrum activity and investigated their inhibitory effects on the CMV RAP-1 and RAP-2 proteins. Our docking results were encouraging with many phytoalexins exhibiting superior and comparable binding energies with the top-scoring standard plant antiviral compound ningnanmycin towards both the target RAP proteins. Among the phytoalexins, the top-scored compounds glyceollidin (-9.8 kcal/mol) and pisatin (-9.7 kcal/mol) for the RAP-1 target had DG values energetically À7.12 kcal/mol and À5.56 kcal/mol more favorable than the standard ningnanmycin. It is worth mentioning that the energy contribution of the residue Lys523 was higher with stronger H-bond interactions which are essential for the RAP-1 target inhibition. Regarding the RAP-2 target, the top-scored compounds Moracin D (-7.8 kcal/mol) and Piceid (-7.3 kcal/mol) had comparable DG values with the standard ningnanmycin with amino acid residues Lys648 and Arg561 as major contributors forming H-bond interactions for stronger inhibition. The compound sanguinarine had lower DG values which may be due to the poor internal steric clashes and low-resolution structure. Overall, from our computational screening results, it is evident that the phytoalexins have the potential to inhibit the RAP-1 and RAP-2 proteins of CMV and could be further exploited as a lead compound for plant antiviral research. Currently, the proposed lead phytoalexin compounds with anti-CMV potential are being evaluated in our laboratory for their practical application to inhibit CMV replication in vitro and in vivo using anti-CMV bioassay. Additional studies on the proliferation inhibition and the ability to induce systemic acquired resistance (SAR) in crop plants would be interesting to understand the role of phytoalexins in plant immunity. Further studies on their mechanism of action, phytotoxicity, and role in CMV resistance should be investigated for their field-level application in sustainable agriculture as plantbased biopesticides.