Eugenol Ester Derivatives: Synthesis, Insecticidal Activity and Computational Studies †

: Specific structural modifications in eugenol molecule can simultaneously improve the biological activity and reduce side effects of the respective analogues. The esterification of eugenol by two different experimental procedures, and subsequently conversion of one of the esters into the corresponding oxirane was carried out. All derivatives obtained were then evaluated for their effect on the viability of Sf9 ( Spodoptera frugiperda ), cells. In addition, a structured-based inverted virtual screening protocol was employed to identify the potential proteins associated to the observed insecticidal activity. The encouraging results obtained allowed to establish a preliminary structure-activity relationship.


Introduction
Due to the exponential increase in world population, it is necessary to ensure agricultural production that meets the actual food requirements. The improvement in the productivity of agricultural crops implies an incessant need to prevent, control and destroy the pests that affect them, achieved through the extensive use of synthetic pesticides. Although synthetic pesticides represent a plausible approach, they present a serious threat because their uncontrolled use causes negative impacts on the environment (pollution and loss of biodiversity) and on human health [1,2].
Natural products are good alternatives, due to the structural diversity and associated biological activity, so as they are a rich source of inspiration in the design and optimization of active principles in the development of formulations, highlighting the crucial role of plant extracts [3,4]. In this category, essential oils fit perfectly, exhibiting a broad-spectrum of actions, including antibacterial, antifungal, insecticidal and antioxidant activities, as for example eugenol [5,6].
Considering these facts, and as a continuation of our recent interests in alternative pesticides, eugenol derivatives were obtained through esterification and epoxidation reactions and evaluated for their effect on the viability of Sf9 cells. A structure-based inverted virtual screening protocol was employed to identify the potential proteins associated to the observed insecticidal activity.

Biological Activity of Compounds 2a-c and 3 in Sf9 Insect Cells
Aiming the evaluation of the insecticidal activity of the synthesized eugenol derivatives 2a-c and 3, studies were carried out in Spodoptera frugiperda (Sf9) cells, a common pest widely used on the screening of insecticides. For benchmarking purposes, the insecticide chlorpyrifos (CHPY) was used at the same concentration (100 µ g/mL). As can be seen in Figure 1, it is clear that the esterification of eugenol with a nitrobenzene group pontentiate eugenol toxicity, derivatives under study displaying equivalent (compound 2a) or even higher (compound 2b and 3) toxicity than the commercial insecticide, CHPY ( Figure 1). Noteworthy, when the nitro group linked to the benzene ring (compound 2b) is replaced by a methoxy group (compound 2c), the cytotoxicity is completely lost. On the other hand, oxirane formation (compound 3) lead to a slight increase in toxicity ( Figure 1). Viability of Sf9 insect cells exposed to the molecules under study 2a-c and 3 (100 µ g/mL), or medium (control) or the reference insecticide chlorpyrifos (CHPY). *** p < 0.001. Table 2 presents the average scores obtained for the four eugenol derivatives for each potential target calculated with each SFs. Regarding the difference in the values, it must be stated that different SFs are based on different scales and metrics. The score for all the GOLD scoring functions is dimensionless with a higher score yielding a better binding affinity. Vina, on the other hand, uses a metric that is a more precise approximation of binding free energy, meaning that a more negative value is equivalent to better affinity. Generally, the results show good consistency between SFs, with odorant binding proteins, acetylcholinesterases, chitinases and beta-N-acetyl-D-hexosaminidase yielding better scores. On the other hand, targets such as octopamine receptor, N-acetylglucosamine-1-phosphate uridyltransferase (GlmU) and voltage-gated sodium channels, consistently present lower scores for across all the SFs.

Inverted Virtual Screening Results
From each set of targets, the structure with the best score was selected and ranked from the best target to worst, according to the predictions of the different SFs. The overall ranking is listed in Table 2. Globally, considering the results obtained with the several SFs, odorant binding proteins are the most likely target with the highest affinity towards eugenol derivatives, followed closely by acetylcholinesterases. The discrepancy in some of the values of the different SFs, can be explained by the own nature of each SF, as they consider different aspects of protein-ligand binding.
The hypothesis formed is that eugenol and eugenol derivatives can be used as repellents because they can bind to odorant binding proteins or used as pesticides, inhibiting insect acetylcholinesterase.
Interestingly, in the PDB database there is a structure of an odorant binding protein bound to eugenol Apis mellifera (PDB: 3S0E) [7]. This might be an important indicator of the increased affinity of eugenol derivatives against OBPs.

Molecular Dynamics Simulations and Free Energy Calculations Results
In order to validate the inverted VS predictions, molecular dynamics simulations were then performed for the eugenol derivatives complexes formed with the two groups of targets predicted at the inverted VS stage: odorant binding proteins and acetylcholinesterases. The structure with the best score from each group was selected (3K1E for OBP and 1QON for acetylcholinesterases-AChE). The results are detailed in Table 3. The OBP-eugenol derivatives complexes are very stable throughout the simulation and presented an average protein RMSD around 2 Å . The prediction from the inverted VS were confirmed as the ligand RMSD is very low. For AChE-eugenol derivatives, however, the average RMSD is higher, indicating that the system shifted to a more stable conformation in the beginning of the simulation. Also, the inverted VS predictions were validated for this target, as the average ligand RMSD values are bellow or equal to 1 Å .
The average SASA and percentage of potential ligand SASA buried indicate the ligand exposure to solvent, so, and increased SASA and a lower percentage of ligand buried, the more solvent exposure. Compounds 2c and 3 are the ones that are less exposed to the solvent and more buried in the binding pocket of OBP. Regarding AChE, the compounds that are less exposed and more buried in the active site are 2a and 3.
Generally, the Gibbs free energy of association was better for OBP-eugenol derivatives than for AChE-eugenol derivatives. Compounds 2a and 3 are the ones that present the strongest affinity toward OBP. Compound 3 is also the compound that presents the strongest affinity toward AChE than all the other eugenol derivatives studied.
When bound to OBP the ligands are mainly stabilized by Trp105, Leu67 and Ile78. When bound to AChE, the main interacting residues are Trp81, Tyr69 and Tyr322.

Evaluation of Viability in Sf9 Cells
As a model, the Spodoptera frugiperda Sf9 cell line was used. Cells were purchased from Sigma-Aldrich (St. Louis, MO, USA) and maintained in Grace's insect medium enriched with 10% fetal bovine serum (FBS) and 1% penicillin/streptomycin (Pen-Strep) at 28 °C. Cells were routinely subcultured as a suspension culture and assays conducted in the exponential growth phase.
For the assessment of viability, Sf9 cells were plated at a density of 3.0 × 10 4 cells/well, followed by incubation for 24 h with the various compounds. After this period, a commercial solution of resazurin was added (Thermo Fisher A13261, final concentration: 1:10) and fluorescence was measured 60 min thereafter.

Inverted Virtual Screening Protocol Optimization
Considering the relevance of the target and year of publication, a search on Scopus was performed using the keywords: Virtual Screening (VS) and insecticide target. Seventeen studies were selected, and thirteen targets chosen for the inverted VS assays. The targets identified are listed in Table 1. Each structure was extracted from the PDB database [25] and was prepared for docking using the Autodock Vina plugin for Pymol [26] with the removal of crystallographic waters and the extraction of ligands to separate files. The saved ligands were later used for active site coordinates and as reference for root mean square deviation (RMSD) calculations. In the absence of crystallographic ligands, the active site coordinates were obtained by selecting the most important active site residues. Re-docking was used as a quality measure, to evaluate the ability of the docking software in reproducing the geometry and orientation of the crystallographic pose.
The docking programs/scoring functions (SF) used were AutoDock Vina [27] and GOLD [28] (PLP, ASP, ChemScore, GoldScore). The protocol was optimized for each protein target and each SF, to minimize the RMSD values.
The optimized parameters for each SF consisted of the coordinates for the docking region centre, docking box dimension or radius, exhaustiveness, search efficiency, and number of runs. Once the RMSD values between poses (crystallographic and docked) were satisfactory (bellow 2 Å ), the optimized conditions were used for the subsequent stages. The molecules were prepared for docking using Datawarrior [29] and OpenBabel [30] and were docked into each structure with all the five SF in study. A ranked list was prepared based on the average scores of each target.

Molecular Dynamics Simulations and Free Energy Calculations
Molecular dynamics simulations were performed on the four eugenol derivatives in complex with the two most promising targets identified from the inverted VS study: Odorant Binding Protein 1 (OBP-3K1E) and Acetylcholinesterase (AChE-1QON). The Amber18 software [31] was used throughout.
The complexes were treated with the Leap module of AMBER [32]. The protein targets were treated with the ff14SB force field [33], while the eugenol derivatives were parameterized using ANTECHAMBER, with RESP HF/6-31G(d) charges calculated with Gaussian [23,34] and the General Amber Force Field (GAFF) [35]. The complexes were placed in TIP3P water boxes with a minimum distance of 12 Å between the protein-surface and the side of the box and periodic boundary conditions were applied. Counter-ions (Na + ) were added to neutralize the overall charge and the complete systems.
To remove clashes prior to the MD simulation, four consecutive minimizations stages were performed with maximum of 2500 steps. Subsequently, the minimized systems were then subject to an equilibration procedure, divided into two stages: in the first stage (50 ps), the systems were gradually heated to 298 K using a Langevin thermostat at constant volume (NVT ensemble); in the second stage (50 ps) the density of the systems was further equilibrated at 298 K. Lastly, the production runs were performed for 100 ns, in a NPT ensemble at constant temperature (298 K, Langevin thermostat) and pressure (1 bar, Berendsen barostat). A 10 Å cutoff for nonbonded interactions was used along with the SHAKE algorithm, to constrain all covalent bonds. An integration time of 2.0 fs was applied. The final trajectories were analyzed using the cpptraj tool [36] and VMD [37], to confirm that all the systems were well equilibrated. The last 70 ns of the simulation were considered for hydrogen bonding analysis, and cluster analysis of the conformations generated.
In order to estimate the binding free energies of the protein-eugenol derivatives complexes, the molecular Mechanics / Generalized Born Surface Area method [38] was applied using the MM/PBSA.py [39] script from amber. The salt concentration applied was 0.100 mol dm −3 . From each MD trajectory, a total of 1400 conformations were taken from the last 70 ns and the contribution of the amino acid residues was estimated using the energy decomposition method.

Conclusions
In this work, three esters derived from eugenol and the corresponding oxirane from one of these esters were efficiently prepared. The obtained eugenol derivatives were subjected to biological activity evaluation in Sf9 cell line, in order to predict their potential as natural based insecticides. We identified that the three derivatives esterified with a nitrobenzene were those showing higher potency, in some cases higher than the benchmark used.
In the present study, we report the application of an integrated molecular modelling-inverted virtual screening protocol for a selection of four eugenol derivatives in order to find possible protein targets in which they present insecticidal activity. After the target selection and protocol optimization, the eugenol derivatives were docked into each of the thirteen targets with five different SFs (PLP, ASP, ChemScore, GoldScore, Vina). Eugenol derivates showed an increased binding affinity for odorant binding proteins and acetylcholinesterases. The fact that there is, already, in the PDB database a structure of an OBP bound to eugenol, is a strong suggestion that eugenol derivates, could be used as repellents. Funding: This research was funded by FCT under project PTDC/ASP-AGR/30154/2017 (PO-CI-01-0145-FEDER-030154) of COMPETE 2020, co-financed by FEDER and EU. FCT-Portugal and FEDER-COMPETE/QREN-EU also gave a financial support to the research centres CQ/UM (UIDB/00686/2020), CF-UM-UP (UIDB/04650/2020) and REQUIMTE (UIDB/50006/2020). The NMR spectrometer Bruker Avance III 400 (part of the National NMR Network) was financed by FCT and FEDER.

Conflicts of Interest:
The authors declare no conflict of interest.