Introduction
Begomovirus, the largest genus of single-stranded DNA genome family Geminiviridae, poses a serious threat to worldwide agricultural productivity. The global dissemination of begomoviruses has raised serious concerns due to their devastating effects on various economically important plant families [1]. Begomoviruses are vectored by a single species of whitefly, Bemisia tabaci, and are categorized into two groups based on geographic location: Old World (OW) and New World (NW) [2, 3]. Furthermore, on the basis of genome organization, they are classified into monopartite (having a single genome component) and bipartite (having two genome components) begomoviruses. The majority of OW begomoviruses are monopartite, while NW begomoviruses typically have a bipartite genome [1]. In the OW, just a few bipartite begomoviruses occur, and Tomato leaf curl New Delhi virus (ToLCNDV) is an example. ToLCNDV is believed to have originated from the Indian sub-continent and subsequently spread to different countries in Africa and Europe [4, 5]. Until 2012, ToLCNDV was only found in the Asian countries. A notable occurrence, however, occurred in the western Mediterranean basin in 2014, where ToLCNDV was found infecting zucchini squash, melon, cucumber [6], and courgette in France [7]. Subsequently, it caused severe epidemic outbreaks that had substantial economic losses to cucurbit crops. Furthermore, ToLCNDV was reported in 2015 in Tunisia, infecting cucurbits, and in Italy, infecting zucchini squash. ToLCNDV has been identified in several other countries, including China, Italy, Spain, Laos, Seychelles, Algeria, Greece, France, Nepal, and Türkiye [8-13], and a recombinant variant of ToLCNDV was very recently identified in Italy [14].
Although ToLCNDV was initially identified in tomato plants (Solanum lycopersicum L.) in 1998. Until now, it has significantly increased and has extended its host range to approximately 45 other plant species, including bitter gourd, bottle gourd, sponge gourd, chili, cotton, cucumber, okra, potato, pumpkin, muskmelon, melon, squash, Solanum nigrum, and various species [8, 10, 12-16] belonging to Cucurbitaceae, Malvaceae, Fabaceae, and Euphorbiaceae plant families. The typical symptoms induced by ToLCNDV in its natural host plants include chlorosis, stunted growth, curling of leaves (upward and downward), vein swelling (thickening), and necrosis.
ToLCNDV is a typical bipartite begomovirus comprising DNA-A and DNA-B. DNA-A encodes six to seven ORFs, including two to three on the sense strand: AV1 (coat protein [CP]), AV2 (precoat protein), and AV3 (its function not ascribed yet). Four ORFs are encoded on the complementary strand: AC1 (Replication-associated protein [Rep]), AC2 (transcriptional activator protein [TrAP]), AC3 (Replication enhancer protein [REn]), and AC4. DNA-B, on the other hand, encodes two ORFs: BV1 (nuclear shuttle protein [NSP]) and BC1 (movement protein). Both have inter-twined functions of systemic movement of virus and infection of adjacent cells [17]. The functions of all these ORFs have already been studied and discussed extensively [18, 19].
Among virion-sense ORFs of ToLCNDV, AV1, AV2, and AV3 partially overlap with each other (Fig. 1a)). These ORFs have garnered substantial attention because of their important role in pathogenesis and interaction with the host plant’s defense response. AV1 ORF is more conserved in nature, whereas AV2 and AV3 ORF regions are more prone to recombination, contributing to the diversity among NW and OW begomoviruses [20] and playing a crucial role in the emergence of new strains such as new ToLCNDV-ES and ToLCNDV-Le in the Mediterranean region [14]. AV2 has been found to play a crucial role in symptom recovery of the host plant [21], and the interactome of the precoat protein elucidates the mechanism of disease progression in plants [22- 24]. The AV2 gene knockdown mutant of ToLCNDV showed highly reduced asymptomatic infection in Nicotiana benthamiana plants [19]. The interaction of ToLCNDV-encoded AV2 with RNA-dependent RNA polymerase 1 (RDR1) of N. tabacum plays a critical role in the development of symptoms in infected plants, indicating its pivotal role in symptom development [25]. The RDR family encompasses RDR1, RDR2, and RDR6 proteins. These proteins play a key role in the plant gene silencing mechanism and protect the host plants by triggering the immune response against foreign pathogens. Particularly, RDR1 orchestrates the production of siRNAs that silence the viral genes. This phenomenon furnishes antiviral protection in a wide range of plant species, such as Oryza sativa [26], N. tabacum, Arabidopsis thaliana, and Capsicum annum [27]. The precise function of the AV3 gene remains an enigma, with only a single study exploring its interaction with the Swa5 gene [28]. Furthermore, it has only been demonstrated within the ToLCNDV genome organization without any functional characterization [29, 30]. However, it has been anticipated that it interacts with RDR1, thereby potentially modulating the host defense response and influencing the severity of viral infection.
Understanding the molecular interactions of ToLCNDV-encoded AV2 and AV3 with S. lycopersicum- encoded RDR1 is critical for understanding the mechanisms that drive virus-host interactions and the development of disease symptoms. The goal of this study is to employ computational approaches to decipher the underlying protein-protein interactions of AV2 and AV3 with RDR1. Insights gained from these studies may contribute to the development of strategies for managing and controlling ToLCNDV and other related begomoviruses.
Materials and Methods
Sequence Retrieval
Initially, a total of 635 full-length ToLCNDV DNA-A genome sequences were obtained from NCBI GenBank (https://www.ncbi.nlm.nih.gov/, accessed on 14 September 2022). These sequences were subjected to multiple sequence alignment via MEGA 11 using the MUSCLE algorithm [31]. From the aligned sequences, the ones containing AV2 and those AV3 were identified and, among these, a total of four, two with AV3 (accession No. U15015 and AM947506) and two without it (KT948072 and MH577696) were used for the downstream analysis.
Subsequently, the AV2 and AV3 gene sequences of the selected ToLCNDV were translated through the Expasy server and aligned through ClustalOmega, and the alignment image was color-coded using 50% conservation of amino acids (aa) using JalView (Fig. 1b)). Additionally, the protein sequence of RDR1 (Solyc05g007510.2.1) was retrieved from the Sol Genomics Network Database (https://solgenomics.net/, accessed on 08 November 2022).
3D Structure Prediction and Validation
To determine the tertiary protein structures of AV2 and AV3, the I-TASSER, an online server and powerful computational tool capable of automatically generating accurate 3D structural models and predicting the biological functions of proteins based solely on their aa sequences and scores the model on the basis of C score, (https://zhanggroup.org/I-TASSER/) was used [32]. The tertiary structure of the RDR1 was predicted through Alphafold2 Protein Structure and Prediction Database (https://colab.research.google.com/ github/sokrypton/ColabFold/blob/main/AlphaFold2. ipynb#scrollTo=kOblAo-xetgx). To ensure the reliability of the predicted best model, it was further validated using a Ramachandran plot through the PROCHECK server, accessed from SAVES v6.0 (https://saves.mbi.ucla.edu/). The 3D protein structures for all AV2 accessions exhibited a high degree of structure similarity. Likewise, AV3 accessions showed similar 3D structures. Thus, only two proteins, one for AV2 (NCBI accession no. MH577696) and one for AV3 (NCBI accession No. AM947506), were included in the final analysis.
Fig. 1.
The genome organization of ToLCNDV isolates present in Asian countries is shown. These isolates have an additional AV3 ORF that is not found in isolates from Mediterranean countries a). A pairwise amino acid alignment of the AV2 and AV3 proteins from the isolates used in the study is shown b). The colored bars (yellow to brown) indicate conserved residues at a scale of 50%. Amino acids that are present in both proteins are highlighted and show black peaks in the conservation bar.
Protein-Protein Docking
The protein-protein interaction was estimated using the High Ambiguity Driven Docking (HADDOCK v2.2), an information-driven flexible docking approach, which performs the docking process by the information of interface region of the component collected through NMR experimentation or mutagenesis data [33]. To ensure accurate interactions, the binding site of the interacting proteins was added to the active residues in the HADDOCK server. The residues of AV2 and AV3 were predicted through the COACH ligand binding site prediction server (https:// zhanggroup.org/COACH/). The residues of RDR1 were identified through homologs from Protein Data Bank (PDB id: 7EU0_M), which interact with other proteins. The passive residues were set to default, representing buried surface solvability. The interaction was further analyzed using the Ligplot program [34], and 3D images were generated with PyMOL to visualize the interaction.
Molecular Dynamic Simulations
Molecular dynamics simulation (MDS) was used to assess the stability of protein or protein-ligand complexes under different conditions. Molecular dynamics analysis within the biomolecular system was conducted using GROMACS (ver 5.0.2), a free, highly throughput, and open-source software suite for high-performance MDS analysis. The solvation of both complexes was achieved using the TIP3P model. The topology files of the proteins were created through the CHARMM 36 force field, containing information about bonded and non-bonded atoms. CHARMM is a versatile molecular simulation program offering a comprehensive suite of energy functions, enhanced sampling methods, and multi-scale techniques for studying many-particle systems. The system was neutralized by adding Na+ and Cl− ions, and the temperature and pressure were maintained at 300 K and 1 atm, respectively. Each complex of AV2-RDR1 and AV3-RDR1 was simulated for a time of 100 nanoseconds (ns). Several dynamic parameters were computed to study the behavior of the protein complexes, including binding energy, Root mean square deviation (RMSD), Root mean structure fluctuation (RMSF), Radius of gyration (Rg), and Solvent accessible surface area (SASA) [35]. The contact residues were studied and analyzed using the Ligplot tool, which was integrated with PyMol software.
Binding Energy Calculation
To identify the interaction values of best-docked complexes, the energy (∆G) was calculated using molecular mechanics with generalized born and surface area solvation (MM-GBSA) analysis [36]. The application of MM-GBSA calculations serves to seamlessly combine high-throughput MDS with the determination of binding free energies in protein-ligand interactions. This encompasses a spectrum of interaction-free energies, including van der Waals forces, electrostatic interactions, polar solvation effects, SASA, and overall binding energies. It helps in describing the structural stability of the ligand in the protein. The energy Equation is as follows:
∆G = ∆G_complex − ∆G_ligand − ∆G_receptor
Where ∆G is the binding energy, ∆G_complex is the energy of the protein-ligand complex, ∆G_ligand is the energy of the ligand, and ∆G_receptor is the energy of the receptor.
Results
Data Acquisition, Analysis, and Phylogeny
In this study, we initially retrieved 635 ToLCNDV DNA-A sequences from various countries, including Bangladesh, Cambodia, France, Iran, Indonesia, India, Italy, Morocco, Pakistan, Spain, Taiwan, Thailand, and Tunisia. Through multiple sequence alignment and ORF finder tools (https://www.ncbi.nlm.nih.gov/ orffinder/), we observed that these sequences comprised two major groups. The first and largest group of isolates has an additional ORF, referred to as AV3, which is least defined and characterized. The second and smaller group lacked this ORF. Our preliminary investigation unveiled that 503 of these sequences contained the AV3 ORF, and 132 did not. Strikingly, almost all 503 isolates were reported from Asian countries, barring two exceptions: an isolate from Spain (accession No. MW219700) and another from Seychelles (accession No. MH511991). The presence and absence of AV3 in selected ToLCNDV sequences are presented (Fig. S1). The AV2 ORF starts at nucleotide (nt) position 119 and ends at 457, consisting of 339 nt and encoding 112 amino acids (aa), while AV3 ORF is usually located between nt position around 41 and 397, comprising 357 nt and encoding 128 aa. AV3 ORF was found to be present in the oldest known ToLCNDV sequence (accession No. U15015) and was predominantly present in almost all the DNA-A sequences isolated from the OW, with just a few exceptions. Conversely, all the DNA-A isolates reported from Africa and Mediterranean countries lacked AV3 ORF, with only a few exceptions.
Furthermore, a phylogenetic maximum likelihood (ML) dendrogram was constructed using a set of randomly selected sequences. This set comprised ten sequences with AV3 and ten without, ensuring the inclusion of at least one sequence from both Asian and Mediterranean countries. The phylogeny classified the isolates into two clades (Fig. S2). The first clade included all the Asian isolates with the AV3 ORF, while the second clade comprised European and African isolates lacking the AV3 ORF. Based on these findings, two unique accessions, one of which containing AV3 (accession No. AM947506) and one lacking AV3 (accession No. MH577696), were included in the study.
Protein Sequence Analysis, Structure Prediction, and Validation
The protein sequences of AV2 and AV3 consisted of 112 aa and 128 aa in length, respectively, and exhibited minimal similarity, sharing only 14 common aa (Fig. 1b)). Predictions of the structures for both proteins were conducted using I-TASSER with model selection based on the best C-score and TM score to assess the quality and topology of the predicted models. The optimal model for the AV2 protein had a C-score of ‒2.92 and a TM score of 0.38, with fewer helical and coil structures (Fig. 2a)). Conversely, the optimal model of AV3 displayed a C-score of ‒4.12 and a TM score of 0.28, and this model contained over 50% helical elements and coils (Fig. 2b)). The 3D structure of RDR1 was predicted through Alphafold, and the pLDDT plot showed a high confidence score of the structure (Fig. 2c)).
Validation of all predicted structures was accomplished using the Ramachandran plot (Fig. S3). For both AV2 and AV3 proteins, 96.3% and 95.9% of aa, respectively, fell within permissible regions, underscoring the accuracy of the predicted models. Similarly, in the case of the RDR1 structure, 99.7% of aa were within the acceptable region, further validating the accuracy of the predictive model. The scores of all protein structures are summarized in Table 1.
The COACH analysis revealed that both proteins had the propensity to interact with RDR1 by sharing several common binding sites. Although the residues involved in their interaction varied, there were common residues shared with RDR1, namely E825, D877, K880, K881, and E884. The interacting aa residues of RDR1 were discerned via a homologous analysis from PDB (accession No. 7EU0_M), considering its interaction with another reported protein. The identified residues ranged from 867 to 890. During this analysis, passive residues were maintained with their default settings, representing buried surface solvability.
Fig. 2.
The best predicted model of ToLCNDV-encoded AV2 protein a) and AV3 protein b). The best predicted model of RDR1 protein by Alphafold2, showing its per residue confidence metric (pLDDT) structure c). pLDDT showed a confidence score on the basis of color, and blue color showed a higher confidence score. The secondary structures, α-helices (red), β-sheets (green), and strands (grey) in the predicted RDR1 structure are shown d).
Protein-Protein Docking
Docking was conducted using HADDOCK, and the predicted final models were clustered based on pairwise backbone RMSD at the interface. These clusters were then analyzed in terms of their average interaction energies, including electrostatic energy (Eelec), van der Waals energy (Evdw), desolvation energy (EAIR), buried surface area, and Z score (Table 2).
The alignment of both the docked complexes was superimposed to decipher the positioning of AV2 and AV3 on the RDR1 protein (Fig. 3). Notably, both AV2 and AV3 docked at a similar binding site of the predicted RDR1 structure (Fig. 3, 4a), 5a)). Intriguingly, the AV2-RDR1 complex had a higher Haddock score (‒97.6±3.2) and Eelec (‒377.6±36.7), indicating a more robust binding affinity compared to the AV3-RDR1 complex. The binding aa residues of AV2 that formed hydrogen bonds in the AV2-RDR1 complex were E11, S12, R50, K52, C84, C86, P87, H88, and A91. While residues partook in the hydrophobic bonds were H14, G15, R50, A51, K52, Y54, P83, C89, and P90 (Fig. 4b) and 4c)). For RDR1, the binding residues spanned across T520, S521, Y522, H761, Q822, E825, P827, P828, P830, D872, M875, D877, K880, K881, and E884 (Fig. 4c)).
Conversely, in the AV3-RDR1 docked complex, the binding residues of AV3 involved in the formation of hydrogen bonds were K26, H27, R34, R37, L49, S51, P52, R53, and R55. While H27 and R55 also formed salt bridge interactions. Consequently, these bindings fostered the formation of hydrophobic bonds by residues W23, A33, S36, K38, and I50 (Fig. 5b) and 5c)). Concurrently, the RDR1 residues bore a resemblance to those of the AV2-RDR1 complex, nonetheless, more residues including Q822, M824, E825, P827, P830, I832, D872, M875, D877, P878, K880, K881, E884, S887, I888, G896, V897, and Q904 participated in this complex (Fig. 5c)). Table 3 illustrates the shared amino acid residues with RDR1, with which both AV2 and AV3 interact, as well as the amino acid residues that are specific to AV2 and AV3.
Table 2.
HADDOCK scores of the docked AV2-RDR1 and AV3-RDR1 complexes.
Fig. 3.
Superimposed image of AV2-RDR1 complex and AV3-RDR1 complex. AV2 protein is shown in red, AV3 in blue, and RDR1 in green.
Fig. 4.
Molecular docking of the AV2 and RDR1 proteins forming a complex a), 3D binding model revealing the amino acids residues interacting with each other b), and 2D interaction of the AV2-RDR1 complex revealing interacting residues and types of bonds exist between them c). H-bonds are shown by black dotted lines, and hydrophobic interactions are shown with eyelash structures.
Fig. 5.
Molecular docking of the AV3 and RDR1 proteins forming a complex a), 3D binding model revealing the amino acids residues interacting with each other b), and 2D interaction of the AV3-RDR1 complex revealing interacting residues and types of bonds exist between them c). H-bonds are shown by black dotted lines, salt bridges by red dotted lines, and hydrophobic interactions are shown with eyelash structures.
Molecular Dynamic Simulation
Using GROMACS, MDS was carried out, incorporating various parameters to assess the structural stability and dynamic behavior. This involved RMSD, Rg, SASA, RMSF, and MM-GBSA analysis to evaluate the simulation results comprehensively.
RMSD Analysis
The RMSD analysis provides insight into the uniform structural stability of the protein-protein complex in a dynamic system over the temporal scale. The overall average RMSD values computed for both the AV2-RDR1 and AV3-RDR1 complexes were 0.58 nm and 0.61 nm, respectively (Fig. 6a)).
Notably, the AV2-RDR1 complex exhibited higher RMSD values between 2-12 ns with a peak difference of ~0.2 nm, after which the structure reached equilibrium at 88 ns with minimal deviations. Subsequently, an abrupt decline in RMSD values of the AV2-RDR1 complex was noted from 88-100 ns with a peak difference of ~0.3 nm. Conversely, the AV3-RDR1 complex experienced an increase in RMSD values post 85 ns. Although the average values of the two complexes did not differ substantially, their behavior diverged after 80 ns. This disparity suggested that the AV2-RDR1 has relatively greater stability, as evidenced by lower RMSD values.
SASA Analysis
The SASA analysis was conducted for both complexes to compute polar and non-polar molecular surface area to comprehend their behavior in the solvent and to identify conformational changes in them [37] (Fig. 6b)). Notably, no significant difference was noted in the SASA values of the AV2-RDR1 and the AV3-RDR1 complexes, which measured 591 nm2 and 598 nm2, respectively. Over the 100 ns simulation, the AV2-RDR1 complex showed a lower SASA value from 0 to 30 ns compared to the AV3-RDR1 complex. Subsequently, both the complexes showed comparable values until 80 ns. However, the AV3-RDR1 complex deviated, manifesting a higher SASA value from 80 to 95 ns (Fig. 6b)).
Fig. 6.
The time evolution trajectories of both the AV2-RDR1 and the AV3-RDR1 complexes. a) Root mean square deviation (RMSD), b) Solvent accessible surface area (SASA), c) Radius of gyration (Rg), and d) Root mean structure fluctuation (RMSF) over the 100 ns simulation.
The analysis revealed no significant structural changes in either of the complexes during solvent interaction. Furthermore, both complexes exhibited consistent and stable behavior in the solvent environment, with minimal fluctuation. Nonetheless, the AV3-RDR1 complex exhibited somewhat more subtle behavior.
Rg Analysis
Rg values were inferred to gauge the compactness of both complexes, representing the distance between the centers of receptor and ligand atoms where the highest energy is transferred. A higher Rg value indicates a less compact and less stable complex, whereas a lower Rg value shows a stable and more stable complex with a higher folded structure. Significant fluctuations in the Rg values of the AV2-RDR1 complex were observed (Fig. 6c)), with an average Rg value of 3.14 nm. Conversely, the average Rg value for the AV3-RDR1 complex was 3.25 nm, indicative of its relatively more stable nature.
The Rg value of both complexes showed apparent fluctuations at different time scales. The Rg value of the AV3-RDR1 complex exhibited an increase from ~35 to ~70 ns, while the Rg value of the AV2-RDR1 complex decreased during the same time interval, showing opposing structural conformation of the two complexes. After ~70 ns, the Rg value of the AV3-RDR1 complex decreased slightly, indicating structural instability compared to the AV2-RDR1 complex, where the Rg value decreased significantly (Fig. 6c)). The Rg analysis demonstrated that AV3 enhanced the compactness of the AV3-RDR1 complex, fostering a closed conformation and stability. Conversely, AV2 reduced the compactness of the AV2-RDR1 complex, showing lesser folding with the RDR1 protein within the complex.
RMSF Analysis
The RMSF analysis unveiled that the average RMSF values were 0.29 nm for the AV2-RDR1 complex and 0.23 nm for the AV3-RDR1 complex, highlighting that the aa within the AV2-RDR1 complex had slightly more flexibility than the AV3-RDR1 complex (Fig. 6d)). In both complexes, the RDR1 residues (ranging from 0 to 1100 aa) exhibited lesser fluctuations, indicating higher stability of its aa residues and ultimately contributing to overall stability in the complexes. Interestingly, the AV2 aa showed higher fluctuations, particularly at the C-terminal, demonstrating highly pronounced flexibility in comparison to the AV3 aa, which exhibited lesser fluctuation in their residues. Interestingly, the binding aa residues of AV2 (mentioned in Protein-Protein Docking paragraph) that played a pivotal role in the AV2-RDR1 complex showed lower flexibility. This reduced flexibility may contribute to the enhanced stability of this complex. The RMSF analysis underscores that there is a lower likelihood of detachment between the RDR1 and AV3 in comparison to AV2, which may be due to the lower prevalence of coil and turn structures in the AV3 configuration.
MM-GBSA Analysis
The MM-GBSA analysis was employed to monitor the relative binding energy (∆G) associated with both complexes over 100 ns simulation. The outcomes of these calculations are presented in Table 4.
The free binding energy of the AV3-RDR1 complex was higher than the AV2-RDR1 complex, measuring ‒59.64 Kj mol‒1 (Table 4). The key residues involved in these interactions displayed a distinctive energy profile, and the graphical representation showed the energy hierarchy of the residues from the lowest to the highest values (Fig. 7).
Notably, within the AV2 protein, aa with the highest binding energy were P87, P90, H88, C89, C86, C84, S12, G15, H14, and P83 from a higher to a lower energy profile (Fig. 7a)). These aa played an important role in the interaction of the AV2-RDR1 complex and stabilized it. On the other hand, just a few aa, including R53, S51, W23, H27, P52, R34, K49, and I50 of the AV3 protein, were found to interact with RDR1 but with relatively higher binding energies (Fig. 7b)). The residues of RDR1 also played a role in the interaction, with G884 being common to both AV2 and AV3 complexes, showing its significance in the interaction process.
The outcome of the MM-GBSA analysis highlighted that the AV3-RDR1 complex exhibited higher energy, and nine key aa residues played a critical role in the stabilization of the protein-protein complex. The higher binding energy of the AV3-RDR1 complex underscores a stronger interaction, contributing to the comprehension of their biological significance.
Table 3.
Common and unique amino acid residues of RDR1 interacting with AV2 and AV3.
Table 4.
The free energy (∆G) parameters and the binding score of the AV2-RDR1 and the AV3-RDR1 complexes.
| Complex Name | ∆Evdw | ∆Eelec | ∆Gb | ∆GSA | ∆Gbind |
|---|---|---|---|---|---|
| AV2-RDR1 | ‒105.48 | 912.15 | ‒835.08 | ‒14.23 | ‒42.65 |
| AV3-RDR1 | ‒90.29 | ‒264.6 | 309.03 | ‒13.77 | ‒59.64 |
[i] All energies are measured in kJ.mol‒1. Abbreviations used in the table are ∆Evdw = contribution of van der Waals force; ∆Eelec = electrostatic force calculated by MM force; ∆Gb = solvent free energy calculated through GB force; ∆GSA = solvent accessible surface area value; ∆Gbind = free binding energy.
Discussion
ToLCNDV is a serious looming threat to agricultural crops, causing significant losses in both crop yield and economic values. Its uncontrolled proliferation could reap a perturbation in the food supply chain across Asian, African, and Mediterranean countries. ToLCNDV is mainly divided into two main strains, referred to as the Asian and Mediterranean strains, based on their geographic distribution [4]. Moreover, we observed a discernible difference among all ToLCNDV isolates pertaining to the presence and absence of the relatively least explored AV3 ORF. It is a well-established phenomenon that the interplay between viral proteins and host factors is critical in determining plant virus pathogenicity and spread. This study investigated the interaction of ToLCNDV-encoded AV2 and AV3 with the RDR1, as well as the implications for virus infectivity and its spread across Asia and the Mediterranean countries.
Fig. 7.
The binding energies of the AV2-RDR1 a) and AV3-RDR1 b) complexes. The interacting residues of AV2 (red bar), AV3 (blue bar), and RDR1 (green bar), which play an important role in the interacting complexes, are mentioned at the lowest and largest values.
The presence of AV3 ORF may account for the distinct characteristics of these two ToLCNDV strains, as most Mediterranean strains, barring a few, lack this feature. We hypothesized that AV3 might share a common (Viral Suppressors of RNA silencing) VSR role with AV2 and could potentially serve as a virulence determinant, facilitating the enhanced pathogenicity and extensive dissemination of ToLCNDV. A previous study has shown that the expression of RNA silencing suppressors augments the pathogenicity and systemic movement of ToLCNDV [38]. Additionally, this proposition aligns with the fact that all Mediterranean ToLCNDV isolates demonstrate poor adaptation to tomato plants [39]. Our hypothesis is supported by another study, which examined the infectivity of ToLCNDV isolates from the Indian subcontinent and the Mediterranean in cucumber plants. Their investigation demonstrated that Indian isolate induced more severe infectivity in comparison to the Spanish isolate [40].
In this study, the in silico interactions of ToLCNDV-encoded AV2 and AV3 with RDR1 were investigated. Both viral proteins docked to RDR1 at similar positions suggest a convergence in their interactions with this key host factor. This raises the possibility of a conserved mechanism by which AV2 and AV3 exploit RDR1 to enhance viral infectivity. Nonetheless, the AV3-RDR1 complex exhibited remarkably higher docking attributes and scores, including RMSD, Rg, SASA, RMSF, and MM-GBSA. Additionally, this complex was more compact and stable than the AV2-RDR1 complex. Notably, the interaction between AV3 and RDR1 involved a greater number of aa residues, indicative of a robust interaction. In this case, different binding possibilities can be speculated. For instance, AV2 and AV3 might bind to RDR1 at a site with partial overlap, which could result in a competitive interaction for binding sites. Conversely, their binding could follow a sequential pattern, allowing both proteins to cooperate without direct competition. Nevertheless, a thorough and detailed investigation is needed to shed light on these possibilities, making it an intriguing avenue for future research endeavors. It is noteworthy that RDR1 shared analogous binding sites with its homolog RDR2 during its interaction with the NPRD1 subunit of Pol4 [41]. Moreover, RDR1 has also been found to play a role in the TGS pathway in the case of methylation of ToLCGV, which helps the host plant in recovery [25, 42]. The outcomes indicate that the interaction of both AV2 and AV3 proteins with RDR1 could potentially inhibit the RDR1 ability to amplify the production of siRNAs, a crucial component of antiviral defense. This inhibition, in turn, could increase the pathogenicity of ToLCNDV in tomato plants [40]. This observation provides the impetus for speculation that the presence of AV3 might confer enhanced interaction capabilities with tomato plants, potentially contributing to the widespread prevalence of ToLCNDV in the Asian continent while potentially accounting for its poor adaptation to tomato plants in Mediterranean countries [39]. This key query formed the pivotal rationale for undertaking this study. However, to validate the in silico findings, it is essential to conduct an in vivo protein-protein interaction study using either the yeast-2-hybrid system or immunoprecipitation.
A mutation in the AV2 gene of ToLCNDV resulted in attenuated symptoms in N. benthamiana plants. When the DNA-A harboring the AV2 mutation was inoculated alone, it produced greatly reduced asymptomatic infection. When inoculated with DNA-B, the mutation did not affect symptoms or infectivity but did delay the onset of symptoms [19]. The delayed and attenuated symptoms may be linked to the absence of AV2’s role in gene silencing. This role may be complemented by either AV3, when DNA-A was inoculated alone, or by other suppressors encoded on DNA-B, which efficiently complemented the missing functions of gene silencing and movement of AV2. However, this is far from conclusive and warrants further investigation.
The discovery of the AV3 gene in ToLCNDV has significant implications for our understanding of viral pathogenicity, virus-host interactions, and the development of plant resistance strategies. AV3 likely plays a pivotal role in ToLCNDV’s ability to infect plants, evade defenses, and manipulate host cellular processes. Its identification opens avenues for breeding or genetically engineering virus-resistant crops by targeting AV3. A more comprehensive approach could involve simultaneously targeting both AV2 and AV3 using gene editing tools or biocontrol agents to achieve robust control over ToLCNDV. Furthermore, AV3 could serve as a valuable molecular marker for early disease detection, enabling timely management strategies. Its potential applications extend beyond ToLCNDV, offering insights into related viruses and biotechnological advancements, such as RNAi-mediated resistance. Overall, the identification of AV3 provides valuable insights for combating ToLCNDV and other plant viruses.
Here, we speculate the possibility that the encoding of AV3 by ToLCNDV might be linked to a least recognized concept of “dual deception” in viruses and bacteria [43, 44]. Dual deception is a strategy where the virus employs two distinct mechanisms, often through different proteins or pathways, to evade or manipulate the host’s immune system and enhance viral replication. The mechanism of “dual deception” may be well-established in viruses but lacks widespread recognition. Viruses achieve dual deception in two ways: by mimicking the host’s own molecules and by inhibiting the host’s immunological response. In the latter mechanism, viruses present two different faces to the immune system: one that is harmless and one that is pathogenic. We hypothesize that AV3 is harmless, unlike AV2, and it is not recognized by the host immune system [45] but performs the same function as AV2. However, this speculation warrants a thorough investigation to draw a concrete conclusion.
Conclusions
Conclusively, this study highlighted the complex interplay between ToLCNDV-encoded AV2 and AV3 proteins and tomato RDR1. The findings imply that despite minimal amino acid sequence homology and distinct 3D structures, both AV2 and AV3 differentially bind to a common binding site of RDR1, suggesting a shared pathway to suppress host plant defense to augment the viral infectivity. Additionally, the study suggests that AV3 establishes a strong interaction with tomato RDR1, offering speculative insights into why Asian isolates of ToLCNDV are better adapted to tomato plants than Mediterranean isolates. Exploring AV3 can help researchers gain a better understanding of ToLCNDV pathogenicity, replication, movement, and host interactions. Such insights hold promise for devising efficacious strategies in disease management and crop protection strategies.
