Potential Acetylcholinesterase Inhibitor Acting on the Pesticide Resistant and Susceptible Cotton Pests.

Gossypium spp., produces economically important cotton fiber, and its yield is highly affected due to pest attacks. Insecticidal target site mutation is one of the reasons behind insecticide resistance to a wide range of pesticides. Acetylcholinesterase (AChE) protein sequences from major pests of cotton were analyzed to assess various physicochemical properties, presence of motifs, and understand evolutionary relationship. The impact of three mutant AChE1, A. lucorum A216S, B. tabaci F392W, and A. gossypii A302S, on the strucutral stability was assessed, and F392W_AChE1 was selected based on 100 ns molecular dynamics simulation. Virtual screening of the zinc database and high-throughput virtual screening, standard precision, and extra precision docking resulted in the identification of six compounds. The six identified compounds and six known commercial pesticdes were docked with three mutant and three wild type AChE1, and one (C1) was selected based on Tice criteria. The conformational and interaction stability of the AChE1-C1 and F392W_AChE1-C1 complexes were monitored at 100 ns Gromacs simulation and were found to be thermodynamically favorable. Therefore, C1 may have the potential to bind to the resistant and susceptible strains of cotton pest, and the resistance developed by insects could be arrested. Furthermore, synthesis and field study of C1 will lead us to a better understanding of the efficacy of the identified compound.

Introduction

Sustaining food security on par with the ever-growing population is a challenge to agricultural productivity. Modern farming practices have resulted in a better yield from various crops of social and economic importance.[1] However, there are still challenges such as dependence on rainfall, pest attack, and soil erosion, among others.[2]

Cotton fiber producing the Gossypium species holds commercial significance globally. India is the largest producer of cotton (29 million bales, 2019–2020), but its productivity/hectare achieved is not on par with other nations.[3] One of the prominent reasons for the reduced yield of cotton is insect pests and crop diseases. Effective control of insect pests is crucial for better crop yield and also in containing insect-borne conditions.[4] Pesticides such as organophosphates (OP) and organochlorines (OC) are widely used to control insect pests by seed treatments and foliar applications.[5] Prolonged and intensive pesticide application has resulted in pest resurgence and pesticide resistance.[6] Over 500 studies have reported a decrease in the effect of insecticides due to insecticide resistance.[7]

Genetic modifications such as transcriptional changes and point mutations in coding regions render insecticide resistance to insects. These allow higher rates of insecticide detoxification and result in target site insensitivity.[8] Another important mechanism for insecticide resistance is metabolic resistance. Detoxification enzymes in insects can block or hydrolyze the insecticide before it reaches the target site.[9]

Acetylcholinesterases (AChE) are catalytically active enzymes belonging to the multigene carboxylesterase family found across the species.[10] Insecticides such as OP and carbamates usually phosphorylate or carbamylate the serine residue present in the active pocket of AChE to inhibit its activity.[11] This process results in the accumulation of acetylcholine in the synapses, leaving the acetylcholine receptor permanently open, increasing the excitement in nerves, leading to insect death.[12] An altered AChE insensitive to the active metabolite can escape from the adverse effects of insecticides.[13]

Extensive loss of biodiversity and land degradation are due to the usage of synthetic pesticides and herbicides. Implementation of organic farming practices is highly recommended as a part of UN Sustainable Development Goals 2030 (SDG).[14] Therefore, plant-based organic chemical compounds are a better option for biocontrol measures.[15]

Molecular modeling approaches of bioinformatics have enhanced the understanding of biological systems. They have also helped in the identification of alternative compounds[16,17] to inhibit targets such as AChE.[18] Several studies have employed high throughput virtual screening,[19] QSAR,[20,21] molecular docking,[22−24] and molecular dynamics[25] for identifying potential inhibitors to AChE, while other approaches have also been used to study the mutations in this enzyme. Structural instability due to the breaking of intermolecular hydrogen bond has been attributed to variations in AChE1.[26] Mutations in AChE analyzed from different insect pests are found to alter the interaction between OP and AChE, thereby affecting the phosphorylation of the enzyme.[22] Pesticides such as OP are proven to induce toxicity in humans by interfering with their metabolic pathways.[27]

Due to the above reasons, it becomes a necessity to find solutions to address the issue of pesticide resistance. Considering the importance of active site-based mutation, we performed a comparative analysis of AChE coding genes from eight major cotton pests, and preference was given to AChE from insecticide resistance strains. Natural compounds having an inhibitory effect on AChE from insecticide-resistant and insecticide-susceptible strains were identified by virtual screening, molecular docking, and molecular dynamics simulation studies.

Materials and Methods

Data Retrieval

A flowchart illustrating the methodology employed in this study is shown in Figure . Insect specific (keywords- insect/insecta) 1710 AChE sequences (keywords, “acetylcholinesterase”, “acetyl cholinesterase”, “acetylcholine esterase”, and “ACHE or AChE”) were downloaded from the NCBI protein database. Cotton pest-specific sequences were used to perform multiple sequence alignment using the Clustal-W module in MEGA7.[28] Redundant sequences (100% identity) were discarded using Elimdup, and 24 (>250 AA) sequences were considered for further study.

Figure 1

Steps involved in the analysis of AChE mutation and pesticide-like compound designing. Color codes: blue: input dataset; green: first phase involving sequence analysis; yellow: second phase involving analysis of the structural impact of point mutation using molecular modeling studies; light gray: docking preparation; blue-gray: XP docking; dark gray: virtual screening procedure; orange: molecular dynamics and interaction analysis.

Multiple sequence alignment involving 24 unique sequences was used to generate a residue-wise similarity graph using Plotcon (http://www.bioinformatics.nl/cgi-bin/emboss/plotcon). A neighbor-joining method-based phylogenetic tree was constructed with 1000 bootstrap replicates based on the Jones–Taylor–Thornton (JTT) model.

3D Structure Prediction and Model Validation

The experimental structure of any cotton-pest specific AChE was not available. So, 3D models of 24 AChE sequences were generated using the Swiss-model workspace server.[29] Modeled structures were validated using the different programs available at the Structural Analysis and Validation Server (SAVES) and the Ramachandran plot using the RAMPAGE server.[30]

Physicochemical Properties and Mutation Analysis

The physicochemical properties were computed using Protparam,[31] putative protein motifs were identified using Multiple Em for Motif Elicitation (MEME) suite,[32] and amino acid compositions were calculated using our in-house Perl script. Mutations in AChE were identified from the literature, and the impact of point mutation on protein stability was ascertained using I-Mutant Server2.0.[33]

Conformational Stability Using Molecular Dynamics Simulation

Structural flexibility and conformational stability of the highly deleteriously mutated (F392W) and wild type structures were studied using molecular dynamics simulation (MDS) for a period of 100 ns using Gromacs version 4.5.3 installed on the High-Performance Computer Cluster (SRM-HPCC). Protein topology parameters were obtained by implementing GROMOS96 force field. Each system was immersed into an SPC water model in a cubical box with a minimal distance of 10 Å from the protein surface to the edge of the box. A periodic boundary condition was applied such that the number of particles, pressure, and temperature remained constant throughout the simulation. Eleven sodium ions were added to neutralize the system. Unfavorable contacts in each system were energy minimized, and the temperature was maintained by applying the Berendsen algorithm.[34] Final production MD was carried out for 100 ns using the isothermal–isobaric (NPT) ensemble. During the simulation, the temperature was maintained constant at 300 K. To treat long-range Coulombic interactions, the particle mesh Ewald method[35] was used. LINCS algorithm[36] was applied to measure bond lengths between hydrogen bonded atoms, and a time step of 2 fs was allowed. Coulombic interactions were truncated at 0.9 nm, and the van der Waals force was maintained constant at 1.4 nm.[37]

Principal Component Analysis

Advanced analysis of trajectories generated from MDS was done using essential dynamics. From principal component analysis (PCA), the motion of the protein was determined by extracting the concerted movement of the protein in different frames during the simulation. PCA was done in two steps: (i) constructing a variance/covariance matrix using C-alpha atoms and (ii) diagonalization of the covariance matrix. PCA was performed using g_anaeig and g_covar tool of Gromacs.[38]

Ligand and Protein Preparation

Zinc natural database (156,601 compounds)[39] and six known cotton pest pesticides (chlorpyrifos (OP), malathion (OP), acephate (OP), methomyl (C), thiodicarb (C) and deltamethrin (P)) structures were downloaded and prepared using the Ligprep module of Maestro 10.2 (Maestro v 10.2 Schrodinger LLC, New York, NY). During ligand preparation, the pH was set to 7.0, and the tautomeric and ionization states were prepared using the Epic module. A maximum of 32 conformers were generated for each ligand structure.[37] The native and mutant AChE proteins were prepared using the protein preparation wizard of Maestro, and the structures were energy minimized using the Optimized Potentials for Liquid Simulations 2005 (OPLS2005) force field. A cubic grid box (gridbox X, Y, and Z ranges were 30, 30, and 30, respectively) was created around the active site comprising SER262, HIS501, GLU388, and TRP392 of AChE protein.

Virtual Screening

Mutant protein (Genbank acc: ABV45412.1, Bemisia tabaci F392W mutation) was used for further study. Virtual screening and Tice rule[40] based filtering was employed to identify compounds that could have the potential to inhibit AChE in both insecticide-resistant and insecticide-susceptible strains. The Glide docking program was used to perform virtual screening in the following order: high-throughput virtual screening (HTVS), standard precision (SP), and extra precision (XP) docking. HTVS docking was performed with a flexible docking algorithm against BT_F392W eliminated zinc compounds that had more than 300 atoms and 50 rotatable bonds. The HTVS docked compounds were filtered based on the Tice rule. The top 10% of compounds based on binding energy from HTVS were docked in the SP mode, and the top10% from SP screening were docked in XP mode. Molecules with the best Glide score and Glide energy were visually inspected and considered for further analysis.

Docking

The top six hits obtained from virtual screening and six known pesticides (discussed earlier) were docked with three wild type and their respective mutant variants. This was done to ensure that the identified compounds have affinity to both the insecticide-susceptible and insecticide-resistant AChE (all three reported mutations).

Stability Analysis of Docked Complexes

The protein target structures (mutated and wild type) were simulated as mentioned earlier. The simulation of docked complexes was performed as follows. One compound (identified from virtual screening studies) having a good affinity toward all three known mutations and its respective wild type was taken for further analysis. The complexes involving the identified compound and mutated AChE F392W and wild type were simulated, and the stability was monitored for 100 ns using Gromacs version 4.5.3. The ligand topology was obtained using Prodrg webserver,[41] and other parameters used were the same as used in the protein simulation discussed above.

Result

An initial multiple sequence alignment (MSA) of 54 cotton pest AChE sequences helped us in removing duplicates. Thirty-six unique sequences comprise Aphis gossypii (12), Helicoverpa armigera (7), Apolygus lucorum (6), Lygus Hesperus (3), Bemisia tabaci (3), Spodoptera littoralis (2), Spodoptera exigua (2), and Spodoptera litura (1). These sequences were aligned, and the similarity profile of aligned sequences obtained from Plotcon is shown in Figure . It can be seen that higher similarity region was found around ∼260 to ∼510 residues, and the area ∼300–310 and ∼400–420 seems to be relatively divergent.

Figure 2

Similarity profile of AChE aligned sequences. The horizontal axis represents the relative residue position and the vertical axis represents the similarity.

Phylogenetic Tree

The NJ tree generated by MEGA was analyzed to gather the evolutionary information of the AChE enzyme (Figure ). The tree has two main clades: AChE1 and AChE2. Both the clades display an aggregation of sequences under their respective species-wise sub-clades involving species Aphis gossypii, H. armigera, B. tabaci, and A. lucorum. It could also be understood from the tree that both the AChE1 and AChE2 sequences of L. hesperus accommodate themselves in close relationship with A. lucorum. From the tree, it is noted that AChE sequences of the same order and species are highly conserved. This provides us the base for developing broad-spectrum pesticides having the potential to inhibit AChE in pests from the same order.

Figure 3

Phylogenetic Tree of AChE. Clade 1 contains AChE2 and clade 2 contains AChE1 sequences.

Primary Sequence Analysis

The length of AChE sequences ranges from 300 to 700 aa. The amino acid composition is represented as a chart in Figure SF1. Physiochemical properties are given in Table ST1. Physiochemical properties showed that negatively charged residues (Asp+Glu) were higher than positively charged residues (Arg + Lys). Hydropathicity (GRAVY) shows that the majority of the sequences are hydrophilic. The average molecular weight of AChE is around 60 kDa. All the sequences contained a motif related to the carboxylesterase family (Figure SF2). Motif profiling was done to ensure that the selected sequences are AChE indeed.

Mutation Analysis

The molecular basis of insects showing resistance to organophosphates, pyrethroids, carbamates, and chlorpyrifos is due to mutations in the insecticidal target site. From the literature, three-point mutations in AChE1 were identified in A. lucorum (AOS89450.1, hereafter AL_MT), B. tabaci (ABV45412.1, hereafter BT_MT), and A. gossypii (BAD51410/BAD51412, hereafter AG_MT) at A216S, F392W, and A302S, respectively, and these details are listed in Table .

Table 1

Mutation Reported in AChE of Different Species

protein ID	gene ID	organism	family	mutation	reference
AOS89450.1	KT805420.1	Apolygus lucorum	Miridae	A216S	(48)
ABV45412.1	EF675187	Bemisia tabaci	Aleyrodidae	F392W	(49)
BAD51410/BAD51412	AB180401, AB180403	Aphis gossypii	Aphididae	A302S	(47)

Stability Analysis of Mutant Proteins

The I-mutant2.0 server predicted the impact of these mutations on enzyme stability. I-mutant calculates the energy difference between native and mutant protein based on Gibbs’s free energy.[42] The difference between the unfolding Gibbs free energy value of the mutated protein and wild type (kcal/mol) gives us the DDG value, which measures the effect of single point mutation on protein stability. Protein with a DDG value less than zero will have reduced stability and vice versa. The DDG of AL_MT, BT_MT, and AG_MT is −0.44, −1.37, and −1.13, respectively. BT_MT F392W had a DDG value that was less compared to the other two mutations. Figure shows mutations occurring in sites in different AChE1. Alanine and serine are categorized as less hydrophobic and hydrophilic amino acids, respectively. Figure A,B shows that A216S and A302S are on the coil region. Phenylalanine is hydrophobic, and tryptophan is relatively less hydrophobic; this mutation lies in the helix of AChE (Figure C).

Figure 4

(A) Aphis gossypii (A216S), (B) Apolygus lucorum, and (C) Bemisia tabaci (A302S). Substitution of A to S replaces one of the methylene hydrogens with a hydroxyl group. F and W are derivatives of alanine, F has phenyl substitution, and W has indole substitution in the beta carbon of alanine. Structural differences due to mutation can be noticed.

Molecular Modeling and Validation

The modeled AChE structures were validated by using ERRAT, Verify3D, and Ramachandran (RC) plots. Details of the templates and validation report are shown in Table . The ERRAT score of all structures was greater than 80, and in the RC plot, around 98% of the residues were found to be aggregating in the allowed region. This validation confirms that the modeled structures were reliable for further study.

Table 2

Homology Modeling: Template Details and Validation Report

		template details					RC plot
query	query organism	PDB ID	query coverage	identity %	verify 3D	ERRATscore	% favored	% allowed	% outlier
AOS89450.1_MT	Apolygus lucorum	2W6C_X	90	45	90.96	92.263	94.5	5.1	0.4
AOS89450.1_WT	Apolygus lucorum	2W6C_X	90	45	90.96	92.263	94.5	5.1	0.4
ABV45412.1_MT	Bemisia tabaci	2W6C_X	82	44	95.51	93.269	94.4	4.9	0.7
ABV45413.1_WT	Bemisia tabaci	2W6C_X	82	44	94.57	94.615	94.7	4.3	0.9
BAD51412_MT	Aphis gossypii	1DX4_A	81	55	93.98	86.667	94	4.9	1.1
BAD51412_WT	Aphis gossypii	1DX4_A	81	55	93.98	86.667	94	4.9	1.1

AChE Structural Analysis

The modeled structure of AChE was analyzed with reference to AChE of Torpedo californica (PDB ID: 2ACE)[43] to gain knowledge about active sites. The active site of AChE contains two subsites, the (i) esteric site and (ii) anionic site. In the esteric site, the Ser–His–Glu catalytic triad is present.[43,44] This catalytic triad is present at sites 262, 501, and 388 in BT_WT/MT, 276, 521, and 405 in AG_WT/MT and 215, 455, and 341 in AL_WT/MT.

The anionic site is formed by four residues Trp–Tyr–Tyr–Phe.[44] The location of the four-residue anionic site in BT_WT/MT, AG_WT/MT, and AL_WT/MT is 147, 193, 391, and 392; 145, 201, 408, and 409, and 100, 146, 344, and 345, respectively. Other sites present in AChE are the acyl binding pocket, oxyanion hole, and peripheral anionic site (PAS). In BT_MT, Phe at 392 of the anionic site is mutated to Trp. In AL_MT, the Ala residue at the oxyanion hole is mutated to Ser. Cartoon representation of AChE active sites is shown in Figure .

Figure 5

Cartoon Representation of the AChE active pocket.

Molecular Dynamics Based Conformational Flexibility and Stability Analysis

I-Mutant2.0 server predicts the change in free energy change. DDG greater than or equal to zero indicates the higher stability of the mutant protein.[33] The F392W mutation has higher negative DGG than others and, therefore, is known to decrease the stability of the protein. So, the structural and functional behavior of wild and mutant (F392W) AChE was studied using the Gromacs MDS package for a period of 100 ns.

Potential Energy

The potential energy observed during the simulation plotted against time is shown in Figure . From the figure, it can be seen that both the systems have stable energy throughout the simulation. The average potential energy reported by wild type protein was −1,316,214 KJ/mol, and the mutant was −1,310,308 KJ/mol. The overall difference in energy between the two types is around 600 KJ/mol. As stated earlier, the simulation confirms that both the systems have attained equilibrium and, therefore, can act as a reference for our further similar study involving docked complexes.

Figure 6

Potential energy of the systems plotted against time. Black and red lines indicate the potential energy of the AChE wild type and AChE mutant type.

RMSD of Protein

Global changes in the protein structure upon mutation were monitored as a function of RMSD plotted against time. The backbone RMSD profile is shown in Figure . The wild type AChE and mutant AChE-F392W structure reached equilibrium after 30 ns. Approximately, after 70 ns, the systems were more converged, and very minor fluctuations were noted in mutant protein at 85 to 90 ns. This minor deviation reflects the changes in protein stability upon amino acid substitution. Also evident from the figure is that RMSD was higher in AChE than AChE-F392W. The maximum RMSD was around 0.37 and 0.35 nm for wild and mutant types, respectively.

Figure 7

Backbone RMSD of the protein structures. The X-axis is RMSD (nm), and the Y-axis is time (ns). Black and red lines indicate the RMSD of the AChE wild type and AChE mutant type, respectively.

C-Alpha RMSF and Radius of Gyration

The wild and mutant AChE showed different fluctuation patterns, and the mutation induced changes in the flexibility of the protein, indicating that the mutation may affect the function of the protein. The RMSF fluctuation centered around the mutated residue, as shown in Figure , emphasizes a similar trajectory in both the mutant and wild type AChE. The radius of gyration (R) chart of wild and mutant AChE against time is shown in Figure . It can be seen that the wild type and mutant AChE were compact throughout the simulation. The average R values of wild type and mutant AChE were 2.27 and 2.43 nm, respectively. This also suggests that wild type AChE is more compact, and the increased R in mutant type AChE could loosen the active pocket, which makes the inhibitor molecule bind on the surface instead of deep inside the binding pocket, leading to unstable interaction as shown in Figure .

Figure 8

RMSF profile of around the mutant site F392W.

Figure 9

Radius of gyration of the protein structures. The X-axis is RMSD (nm), and the Y-axis is time (ns). Black and red lines indicate the RMSD of AChE wild type and AChE mutant type, respectively.

Figure 11

Surface model of AChE wild and mutant structures showing change in the binding pocket. Red: AChE; yellow: pesticide. In the wild type, the inhibitor molecule is bound deep inside the groove, and in the mutant type, the structure of groove was altered and inhibitor lies on the surface of the molecule.

Hydrogen Bonding

One of the main factors responsible for maintaining protein structure stability is hydrogen bonds.[42] Analyzing the number of hydrogen bonds in wild type and mutant AChE proteins is essential to understand the stability between the residues in these proteins. The intramolecular hydrogen bond is shown in Figure . During the final 10 ns, wild and mutant AChE formed a maximum of 370–432 and 364–434 hydrogen bonds, respectively. The average hydrogen bonds observed were 400.63 and 398.86, respectively. Hydrogen bond counts varied in the mutant compared with the wild type AChE. The change in hydrogen bond counts depicts the power of deleterious amino acid in H-bond formation.

Figure 10

Intramolecular hydrogen profile of the wild and mutant AChE proteins. The X-axis is time in ps, and the Y-axis is the hydrogen bond count. Black lines: AChE wild type; red lines: AChE mutant type.

Essential dynamics was used for a border view of dynamic properties with respect to MDS results. The projection of the first two eigenvectors for wild type and mutant AChE is shown in Figure . The covariance matrices of C-alpha atoms for wild type and mutant AChE were 14.78 and 9.97 nm2 respectively. The covariance matrix of C-alpha atoms for wild AChE was higher compared to the mutant. The projection of the first two eigenvectors in phase space for wild type and mutant AChE explained that the cluster of the wild type is stable, and in mutant AChE, it is expanded in the conformational space due to flexibility.

Figure 12

Principal component analysis. Black and red lines indicate AChE wild type and AChE mutant type, respectively.

Pesticide-like compounds having a binding affinity towards wild and mutant AChE were predicted using the computer-aided pesticide designing approach. The grid of mutant BT_W392F was used for virtual screening against the ZINC natural database in hierarchical mode (HTVS → SP → XP).[45] Tice rule and a three-tier virtual screening protocol (HTVS, SP, and XP) resulted in 39,086 compounds with a binding energy ranging from 4.45 to −10.18 kcal/mol. The top 10 % of HVTS compounds were docked in the SP mode. The SP docked compounds have a binding energy between −1.77 and −10.59 kcal/mol, and finally, 391 SP docked compounds were passed to XP docking, and XP docking resulted in a glide score ranging between −8.46 and −13.11 kcal/mol.

The selected hits from the virtual screening study were re-docked to reaffirm the binding affinity of the compounds toward BT_F392W AChE. The binding energy between AChE and inhibitor molecules and residues involved in forming hydrogen bonds is shown in Table . Residues Tyr184, Gly182, Thr185, Ser262, Asn148, and Glu261 are engaged in H-bond interaction. Trp147 and Trp392 were involved in pi–pi stacking. The compound ZINC95099639 has less binding energy with two hydrogen bonds, and the compound ZINC03812983 has a maximum of three hydrogen bonds with a binding energy of −11.98 kcal/mol. ZINC70455604 has a binding energy of −12.48 kcal/mol with no hydrogen bond interaction with AChE. Figure shows the binding pocket residues of AChE and the residues involved in the interaction with ligand molecules.

Table 3

Binding Energy Details of AChE with Identified Compounds

compound	dock score	energy (kcal/mol)	H-bond	pi–pi stacking
ZINC95099639	–13.11	–31.96	Tyr184, Gly180	Trp147
ZINC03812983	–11.98	–33.84	Ser262, Gly182, Glu261	Trp392
ZINC95099641	–12.11	–35.56	Asn148	Trp147
ZINC70455604	–12.49	–37.56	No	No
ZINC32273181	–12.27	–34.97	Thr185	Trp147
ZINC95098859	–12.12	–32.42	Gly180, Tyr391	Trp147

Figure 13

Binding mode interaction of BT_MT with identified compounds. Green: hydrophobic; purple: charged (positive); cyan: polar; yellow: glycine; orange: hydration set; green line: pi–pi stacking; arrow (purple): backbone H-bond; dotted arrow (purple): side chain H-bond.

Docking with Known and Identified Pesticides

Six known pesticides and six identified compounds were docked to the wild and mutant structures of AL_MT, BT_MT, and AG_MT. The binding energy between the compounds (known pesticide and identified) and AChE (wild type and mutant) is shown in Table . Details of residues involved in hydrogen bond interactions are shown in Table ST2. We found that compound ZINC03812983 had an excellent affinity toward the resistant and susceptible AChE followed by ZINC95098859, ZINC95099641, and others. Only ZINC03812983 (C1) binds to the serine residue present in the catalytic triad of AChE, which is a critical residue in the inhibition of AChE by phosphorylating or carbamylating it. This compound was taken for molecular dynamics study. The interaction of C1 with wild type and susceptible AChE is shown in Figure . The Tice rule property of the identified compounds is shown in Table .

Table 4

Binding Energy Interaction between Native and Mutant ACHE1 and Compoundsa

protein/compounds	ABV45412 (MT) (kcal/mol)	ABV45413 (WT) (kcal/mol)	AOS89450 (MT) (kcal/mol)	AOS89450 (WT) (kcal/mol)	BAD51412 (MT) (kcal/mol)	BAD51412 (WT) (kcal/mol)
chlorpyrifos	–4.63	–4.49	–5.45	–6.36	–5.56	–6.43
malathion	–4.39	–4.61	–4.96	–4.78	–5.45	–6.67
acephate	–3.02	–3.33	–3.46	–3.81	–4.03	–4.28
methomyl	–2.50	–2.92	–3.00	–3.10	–3.13	–3.29
thiodicarb	–1.34	–3.21	–3.21	–3.60	–4.00	–3.06
deltamethrin	–5.47	–6.84	–7.72	–6.86	–10.27	–9.39
ZINC95099639	–13.11	–8.06	NA	NA	NA	–5.71
ZINC03812983a	–12.50	–11.23	–9.91	–9.93	–9.85	–11.69
ZINC70455604	–12.49	–8.35	–9.03	–8.70	–5.86	–7.37
ZINC95098859	–12.41	–9.51	–8.11	–9.58	–8.60	–11.59
ZINC32273181	–12.28	–9.06	–6.95	–10.21	–10.36	–9.18
ZINC95099641	–12.11	–7.68	–11.10	–10.18	–9.29	–14.73

*, best compound; NA, no interaction observed; MT, mutant type; WT, wild type.

Figure 14

Binding mode interaction of C1 with wild type and mutant protein. Green: hydrophobic; purple: charged (positive); cyan: polar; yellow: glycine; orange: hydration set; green line: pi–pi stacking; arrow (purple): H-bond backbone; dotted arrow (purple): side chain H-bond.

Table 5

Tice Pesticide-like Properties of the Compoundsa

compound	name	xlogP	H-bond donors (HBD)	H-bond acceptors (HBA)	molecular weight (MW)	rotatable bonds (RB)
ZINC95099639	4S-methoxylastourvilline	2.25	2	6	357.40	3
ZINC03812983	(−)-heroin hydrochloride	1.1	3	4	286.35	0
ZINC70455604	hydroxy-dimethoxy-methyl-BLAHone	0.57	2	5	328.38	2
ZINC95099641	10-O-demethylcassythicine	4.26	3	5	312.34	0
ZINC95098859	morphine	1.1	3	4	286.35	0
ZINC32273181	mecambroline	3	2	4	296.34	0

Tice pesticide-likeliness property: molecular weight = 150–500, hydrogen bond donor = ≤2, hydrogen bond acceptor = 1–8, rotatable bonds = <12.40

Molecular Dynamics Simulation of Docked Complex

The ZINC compound ZINC03812983 (hereafter referred to as C1) is predicted to have a good affinity toward AChE of both the insecticide-susceptible and insecticide-resistant insects. The stability of the docked complex of wild AChE-C1 and mutant AChE F392W-C1 was studied for 100 ns using Gromacs and molecular dynamics simulation package. The potential energy of both the systems was stable throughout the simulation (Figure ). The highest and lowest potential energies observed over 100 ns for AChE-C1 were −1,339,649 and 1,329,997 kJ/mol, and those for AChE-F392W-C1 were −1,338,953 and – 1,327,752 kJ/mol, respectively. The average potential energies for AChE-C1 and AChE-F392W-C1 were −1,334,895 and −1,333,687 kJ/mol. The average RMSD values of AChE-C1 and AChE-F392W-C1 were 0.25 and 0.39 nm, respectively (Figure ). AChE-C1 has low RMSD, this system was well equilibrated after 70 ns, and the AChE-F392W-C1 has minor insignificant fluctuations after 80th ns.

Figure 15

Potential energy of the systems plotted against time. Black, red, blue, and green lines indicate the potential energy of AChE wild type, AChE mutant type, AChE wild type bound with identified inhibitor C1, and AChE mutant type bound with identified inhibitor C1, respectively.

Figure 16

Backbone RMSD of the protein structures. The X-axis is RMSD (nm), and the Y-axis is time (ns). Black, red, blue, and green lines indicate the RMSD of AChE wild type, AChE mutant type, AChE wild type bound with identified inhibitor C1, and AChE mutant type bound with identified inhibitor C1, respectively.

The RMSF fluctuations of residues surrounding mutation F392W is shown in Figure . From the plot, it can be seen that native AChE has a higher RMSF profile compared to other systems. The RMSF profile of the system changes upon C1 binding, which may result in the loss of function. The radius of gyration is a way to measure the compactness of the system throughout the simulation. The graph shown in Figure shows the compactness of the protein analyzed from the 100 ns trajectory. The average R values observed throughout the simulation were 2.26 and 2.28 nm. A maximum of four and seven hydrogen bonds were observed between AChE and C1 and AChE F392W and C1 throughout the simulation, respectively. Molecular dynamics analysis indicates the strong interaction between AChE and the compound C1.

Figure 17

Central alpha-carbon RMSF of the wild type and mutant AChE proteins. The X-axis is the residue position, and the Y-axis is RMSF (nm). Black, red, blue, and green lines indicate the RMSF of AChE wild type, AChE mutant type, AChE wild type bound with identified inhibitor C1, and AChE mutant type bound with identified inhibitor C1, respectively.

Figure 18

Radius of gyration (R) of the wild and mutant AChE proteins. The X-axis is time in (ps), and the Y-axis is R (nm). Black, red, blue, and green lines indicate the R of AChE wild type, AChE mutant type, AChE wild type bound with identified inhibitor C1, and AChE mutant type bound with identified inhibitor C1, respectively.

Solvent Accessible Surface Area

The solvation effect is another crucial factor that helps to maintain protein stability and folding.[42] Changes in the solvent-accessible surface area of wild and mutant AChE with respect to time are shown in Figure . The average SASA values exhibited by wild and mutant proteins were 224.40 and 219.60 nm2, respectively. The mutant protein showed a relatively lower SASA than the wild type, which indicates the repositioning of amino acids from buried to the accessible area or vice versa.

Figure 19

Solvent accessible surface area. The X-axis is time in ps, and the Y-axis is area in nm2. Black, red, blue, and green lines indicate the SASA of AChE wild type, mutant type, wild type bound with identified inhibitor C1, and mutant type bound with identified inhibitor C1, respectively.

Discussion

Acetylcholinesterase (EC: 3.1.1.7) coded by Ace gene is a crucial enzyme that regulates the neurotransmitter acetylcholine and terminates nerve impulses. Acetylchonline (Ach) is an important neurotransmitter that gets hydrolyzed by AChE at cholinergic synapses, making AChE critical for the normal functioning of the central nervous system.[46] Two main classes of pesticide, namely, organophosphates and carbamate, are known to target AChE.[47]

In the current study, existing AChE sequences of eight cotton pests, which include four species from Hemiptera and four from Lepidoptera, were studied. Phylogenetic analysis of these sequences revealed that the AChE sequence maintained similarity across order and species. Since the AChE sequence is conserved across order (taxonomy level), designing a new class of pesticide (AChE inhibitor) will help control cotton pests of the order Hemiptera. The average sequence length of AChE is around 542 amino acids. Literature survey-based mutation analysis led us to understand three point mutations Bemisia tabaci F392W, Aphis gossypii A302S, and Apolygus lucorum A216S. These three mutations resulted in the shifting of amino acids from hydrophobic to the hydrophilic core, and this might have an impact in the solvent-accessible surface area. Structural analysis of AChE revealed two important subsites in the AChE catalytic site: (i) the esteric site possesses the active OH of serine, imidazole of histidine, and the COO– group of glutamate.[46] These three residues (Ser–His–Glu) form a catalytic triad in which the organophosphates are known to react with serine.[43] (ii) The anionic site is a choline-binding pocket, and it is believed to be a binding site for competitive inhibitors containing the free carboxyl (Asp/Glu) group,[44,46] TRP, TYR, TYR, and PHE. Choline moiety binds to the anionic site mainly through pi–cation interaction.[44]

Conformational flexibility of the wild type and mutant AChE was studied using Gromacs. RMSD, RMSF, hydrogen bonding, R, and SASA trajectories were plotted against time. Both systems have stable potential energy throughout the simulation. The C-alpha amino acid fluctuation reflects the changes in the flexibility of AChE upon mutation (F392W). Principal component analysis of the wild type and mutant F392W AChE also portrayed the impact of this mutation.

Computer-aided pesticide design was carried out to identify compounds having the ability to inhibit AChE of insecticide-resistant and insecticide-susceptible pests. Virtual screening exercise resulted in compounds having an affinity toward AChE. Six compounds identified from the virtual screening study along with six known pesticides were considered for checking their affinity toward AChE from insecticide-resistant and insecticide-susceptible insects. Based on docking studies, ZINC03812983 (C1) is found to have excellent affinity compared to available pesticides. Interaction analysis of C1with mutant and wild type AChE shows that the identified compound C1 binds at the active site of AChE and the OH at the benzene ring forms a hydrogen bond with the Ser262 residue of the catalytic triad. The compound was found to bind in the deep groove of AChE. This compound also forms pi–cation interaction with Trp392, which was a mutant residue (BT_F392W). C1 also binds well with other mutant and wild type proteins and forms a hydrogen bond and pi–cation interaction with serine of the catalytic triad and phenylalanine, respectively. This confirms that the identified compound can bind to AChE of insecticide-susceptible and insecticide-resistant insects. The identified compounds can be taken for further chemical synthesis leading to advanced in vitro and in vivo assessment and product development.

Conclusions

Economically important cotton fiber produced by Gossypium spp. is affected by thousands of pests, which have evolved pesticide resistance to a wide range of pesticides. Mutation in the target’s active site of enzymes such as AChE is attributed to insects developing resistance. In this study, AChE sequences from major pests of cotton were analyzed for its prevalence in these pests and we found AChE1 and AChE2 to be present. Reported mutations in AChE1 were selected, and their impact on the structural and functional properties of AChE1 was studied. We chose F392W_AChE1 and wild type AChE1 for the virtual screening and SP and XP docking of zinc natural compounds. Six compounds from the screening and six known pesticides were further docked with three mutant and three wildt ype AChE1 to select a compound ZINC03812983 (C1) that satisfied the Tice rule and also had binding affinity to the mutant and wild type. A further molecular dynamics simulation of the F392W_AChE1-C1 and wild type AChE1-C1 complexes demonstrated C1 to have firmly bound to both the AChE1 types and displayed acceptable thermodynamic properties and conformational stabililty. Therefore, we believe that C1 has the pesticide-like property and could impact both the insecticide-susceptible and insecticide-resistant species. Further studies such as synthesis and field testing of C1 on different cotton pests are needed and could be taken up as an extension of this study to understand the potential and efficacy.