Unraveling the Mechanisms of Benzimidazole Resistance in Hookworms: A Molecular Docking and Dynamics Study

doi:10.21203/rs.3.rs-4933910/v1

Download PDF

Research Article

Unraveling the Mechanisms of Benzimidazole Resistance in Hookworms: A Molecular Docking and Dynamics Study

https://doi.org/10.21203/rs.3.rs-4933910/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background

Benzimidazole resistance is an emerging challenge among parasitic helminths. It is caused by single nucleotide polymorphisms (SNPs) in specific loci in helminths’ β-tubulin gene. Field studies and laboratory investigations reported resistance-associated SNPs in 4 codon locations with 7 allelic variations among hookworms. This study aimed to determine the effects of these mutations on the binding efficiency and behavior of the β-tubulin protein in four hookworm species against four benzimidazole drugs.

Methods

β-tubulin gene coding sequences of Ancylostoma caninum, A. duodenale, A. ceylanicum, and Necator americanus were retrieved, assessed phylogenetically, and used to construct the 3D structure models of the proteins. The modeled protein structures were verified and edited to contain the reported SNPs: Q134H, F167Y, E198A, E198K, E198V, F200L, and F200Y. Benzimidazole drugs such as albendazole (ABZ), fenbendazole (FBZ), mebendazole (MBZ) and oxfendazole (OBZ) were used as ligands. Molecular docking experiments were performed with the wild-type and mutated proteins. Molecular dynamics simulation assessed the dynamic behavior of the β-tubulin-benzimidazole complex.

Results

In silico docking assessments showed that various amino acid substitutions due to resistance-associated SNPs cause alterations in binding affinities and positions. E198K and Q134H in hookworm β-tubulins substantially weakened the binding affinities and altered the binding positions of benzimidazole drugs. Molecular dynamics analysis revealed that these mutations also caused marked reductions in the binding free energies owing to diminished hydrogen bond contacts with the benzimidazole ligands.

Conclusion

The evidence shown herein indicates that mutations at positions 198 and 134 are detrimental to conferring benzimidazole resistance among hookworms. The presence of these mutations in may alter the efficacy of pharmacological interventions. Hence, further studies should be conducted to assess their emergence among hookworms in endemic areas with histories of chemotherapy.

Parasitology

Antimicrobial Resistance

Deworming

Ancylostoma

Necator

Neglected Tropical Diseases

Hookworm infections continue to pose a significant threat to public and veterinary health worldwide ^1–3. Human hookworm infections are primarily caused by Necator americanus and Ancylostoma duodenale, while A. caninum, A. braziliense, and A. tubaeforme are the primary etiological agents of hookworm infections in animals ^2,3. Ancylostoma ceylanacum is a true zoonotic hookworm widespread in countries within the Asia-Pacific region ⁴. Hookworm infections impose a substantial burden on endemic populations, with an estimated loss of 800,000–4 million disability-adjusted life years ^5,6. Moreover, zoonotic hookworm infections cause approximately 229–500 million cases globally ^7,8. Among animals, hookworms reportedly infect approximately 41% of domestic dogs in the Asian continent, with A. caninum and A. ceylanicum accounting for 27% and 24% of these infections, respectively ⁹. To abate the health threats caused by hookworm infections, numerous control efforts have been mobilized by the public health and veterinary sectors in many countries. Preventive chemotherapy for humans and routine deworming for animals continue to be the cornerstone of hookworm control programs ^10–13. The repeated and sustained administration of benzimidazole drugs, such as albendazole, mebendazole, and fenbendazole, has raised concerns about the development of benzimidazole resistance in soil-transmitted helminths, including hookworms ^6,12–15.

Single nucleotide polymorphisms (SNPs) in the β-tubulin isotype 1 gene result in benzimidazole resistance in helminths ^16,17. Canonical nucleotide mutations F167Y (TTC,TTT to TAC,TAT), E198A (GAG,GAA to GCG,GCA), and F200Y (TTC,TTT to TAC,TAT) cause amino acid substitutions that change the structure of the β-tubulin protein, thereby adversely affecting its binding efficiency with benzimidazole drug ligands ¹⁸. These mutations have been reported in different hookworm species that affect both humans and animals. Recently, multidrug-resistant A. caninum isolates that carry such mutations have emerged in the United States and Canada ^19–21. Moreover, studies on these multidrug-resistant hookworms have resulted in the discovery of novel SNPs that have not been previously found among members of Ancylostomatidae: Q134H (CAA to CAT), E198K (GAG, GAA to AAG, AAA) and E198V (GAG, GAA to GTA ) ^19,22. Likewise, Laboratory investigations that induced resistance in A. ceylanicum have produced inconsistent findings regarding the correlation between these SNPs and benzimidazole resistance ^23,24. One of these laboratory studies identified F200L (TTC,TTT to TTA) mutation which was linked to the induced resistance ²³. Despite these reports, the exact mechanisms underlying benzimidazole resistance in hookworms remain unclear and require further study.

Recent in silico studies determined the key interactions between helminth β-tubulins and benzimidazole drugs. However, these studies focused on the use of albendazole as the sole drug ligand ^25,26. This seems limiting since albendazole is not frequently used against canine or feline hookworms, whereas fenbendazole and oxfendazole are used in deworming preparations ²⁷. Moreover, these studies have not explored the effects of other SNPs reported recently in addition to the canonical mutations. Therefore, to gain novel insights into benzimidazole resistance, this study aimed to determine the effects of all SNPs reported in hookworms on the binding behavior and efficacy of hookworm β-tubulin Isotype 1 protein with several benzimidazole drugs. Here, we report the effects of resistance-associated mutations in the hookworm β-tubulin proteins when docked in silico with benzimidazole ligands.

2.1. Hookworm β-tubulin nucleotide and amino acid sequences

As described previously by Jones et al. ²⁵, four isotypes of β-tubulin exist among the members of Ancylostomatidae (Supplementary Table 1). Among hookworms of public and veterinary health concerns, the coding nucleotide and amino acid sequences of these proteins have only been available for A. caninum, A. ceylanicum, A. duodenale, and N. americanus. These β-tubulin sequences, both nucleotide and amino acids, were obtained from Wormbase ParaSite ²⁸. Moreover, the amino acid sequences were confirmed by searching UniProt ²⁹.

2.2. Phylogenetic Analysis

BioEdit Sequence Alignment Editor version 7.2.5 was used to align both nucleotide and amino acid sequences ³⁰. The sequences were aligned using ClustalW Multiple Alignment ³¹ function in BioEdit. The aligned sequences were trimmed using the alignment editor of MEGA version 11 ³². The best model was determined for both the trimmed nucleotide and amino acid sequences, and phylogenetic trees were constructed. Maximum likelihood trees were drawn based on the K2 (Kimura-2 parameter) + G model for the nucleotide sequences and JTT (Jones-Taylor-Thornton) model for amino acid sequences with 1000 bootstraps in MEGA version 11.

2.3. Homology Modelling and Quality Assurance Checks

β-tubulin proteins of the four hookworm species were modeled using SWISS-MODEL ³³ (https://swissmodel.expasy.org/) [Accessed: 27 May 2024]. The amino acid sequences were used for protein structure modeling, with a previously described β-tubulin (PDB ID no. 6fkj) as the template. 6fkj was selected because it contains β-tubulin proteins (chains B and D), whose crystal structures were determined with a ligand bound to the colchicine binding site in a similar manner to benzimidazole drugs ^34,35. The percentage of sequence identity similarity and GMQE (Global Model Quality Estimate) scores were noted (Supplementary Table 3). Modeled β-tubulins were submitted to various platforms for quality assurance checks. Before downloading the PDB file of the homology-modeled protein, the Ramachandran plot, which gives the percentage of amino acids in their favored conformation in SWISS-MODEL’s structure analysis function, was inspected and noted. Moreover, the modeled proteins were submitted to the UCLA-DOE LAB — SAVES v6.0 (https://saves.mbi.ucla.edu/) [Accessed: 27 May 2024], and the results of the PROCHECK Ramachandran plots and VERIFY3D plots were recorded (Supplementary Table 3).

Since mutations occurring in the β-tubulin isotype 1 gene are the cause of benzimidazole resistance ^16,18, modeled isotype 1 proteins were selected for in silico docking assessments with benzimidazole drugs. Known resistance-associated SNPs at positions 134, 167 198, and 200: Q134H, F167Y, E198A, E198K, E198V, F200Y, and F200L were introduced into amino acid sequences and were used to model proteins using SWISS-MODEL (Supplementary Table 2) ^{19,20,22,23,36}. These mutated proteins were also assessed using PROCHECK and VERIFY 3D (Supplementary Table 4). Both wild-type and mutated hookworm β-tubulins were evaluated for their structural similarity to the mammalian template using the PDB Pairwise Structure Alignment (https://www.rcsb.org/alignment) [Accessed: 29 May 2024]. The mammalian β-tubulin monomer (PDB ID: 1FJJ) was used as the template. RMSD, TM Score, Identity %, Aligned residues, and sequence lengths were documented and assessed.

2.4. In silico Docking

To obtain the best binding poses and their respective affinity values, in silico docking analyses with the modeled β-tubulin isotype 1 proteins and benzimidazole drugs were performed using AutoDock Vina ³⁷ in the PyRx-Virtual Screening Tool version 0.8 (https://pyrx.sourceforge.io/). Before uploading the proteins to PyRx, polar hydrogens were added to all wild-type and mutated β-tubulins using the BIOVIA Discovery Studio Visualizer. These proteins were then used as docking macromolecules for PyRx. 3D ligand structures of Albendazole (as Albendazole sulfoxide, ABZ), Mebendazole (MBZ), Fenbendazole (FBZ), and Oxfendazole (OBZ) were obtained from PubChem ³⁸ (https://pubchem.ncbi.nlm.nih.gov/) [Accessed 25 May 2024]. These structures were imported to PyRx, set as ligands, and energy-minimized.

In silico docking was performed by setting the center of the search grid cube at position 200, with the surface of the cube at 25 Å from the center (coordinates: X = 15.6; Y = 57.5; Z = 33.1). This search grid configuration ensured that all positions where the SNPs occurred were inside the cube. The binding positions with the highest binding affinities (kcal/mol) and the lowest RMSD scores (lower and upper bounds, preferably 0) were considered the best and were saved as PDB files. The binding poses of the in silico docking experiments were visualized using the Receptor-Ligand Interaction function in BIOVIA Discovery Studio Visualizer version 21.1.0.20298 (https://discover.3ds.com/discovery-studio-visualizer-download). Both the 3D and 2D representations of the docking results were evaluated.

2.5. Molecular dynamics simulations

Docked complexes of the top two resistance-associated mutations that caused the highest reduction in binding affinity were analyzed using molecular dynamics (MD) to determine the binding behavior and binding free energies. MD simulations were performed using GROMACS 2023.3. The topology of the benzimidazole ligands (Albendazole and Fenbendazole) was assigned using the Chemistry at Harvard Macromolecular Mechanics (CHARMM) General Force Field (CGenFF) server (https://app.cgenff.com/login). The protein topology of the Ancylostoma ceylanicum β-tubulin receptor was generated using the CHARMM36 force field (July 2022 version), which was incorporated into GROMACS. The protein complex was simulated inside a cubic box with a 10 Å buffer from the box edge to avoid interactions with the periodic image from the adjacent unit cell. The box was filled with TIP3P water molecules and neutralized with Na and Cl ions. Then, 0.15 M NaCl was added. The system was energy minimized using steepest descent for 50,000 steps to remove steric clashes. It was then equilibrated to 300 K and 1 bar using NVT and NPT ensembles, respectively. Both equilibration methods used 50,000 steps or 100 ps of simulation time. Temperature was regulated using the V-rescale thermostat, while pressure was controlled using the Berendsen barostat. The MD run was conducted for 50,000,000 steps or 100 ns, saving snapshots every 5,000 steps, resulting in a trajectory file composed of 1,000 frames. RMSD (Root Mean Square Deviation), RMSF (Root Mean Square Fluctuation), and H-bond plots were constructed to assess the behavior of the β-tubulin-benzimidazole ligand complexes. The absolute binding free energies, determined through MM-GBSA (Molecular Mechanics-Generalized Born Surface Area) and MM-PBSA (Molecular Mechanics-Poisson-Boltzmann Surface Area) methods, were calculated using the entire 100 ns simulations. MM-GBSA and MM-PBSA estimates were derived from the sum of the van der Waals, electrostatic, polar solvation, and nonpolar solvation energies (Eq. 1). Both methods employed the Lennard-Jones and Coulomb potentials for van der Waals and electrostatic energy calculations, respectively. MM-GBSA employed the Generalized Born model to estimate the polar solvation energy, while MM-PBSA used the Poisson-Boltzmann model. Nonpolar solvation energies were assumed to have a linear dependency on the solvent-accessible surface area (SASA).

$$\:\varDelta\:{G}_{bind}=\varDelta\:{E}_{vdw}+\varDelta\:{E}_{elec}+\varDelta\:{G}_{polar}+\varDelta\:{G}_{nonpolar}$$

3.1. Hookworm β-tubulins and their phylogeny

All modeled β-tubulins passed the quality assurance checks based on the constructed Ramachandran Plots in SWISS-MODEL and PROCHECK (Supplementary Table 3). However, the modeled isotype 3 β-tubulins of A. caninum, isotype 4 of N. americanus, and both isotypes 3 and 4 of A. duodenale failed the VERIFY3D check (i.e., the modeled protein scored less than 80%). All isotype 1 β-tubulins of all hookworm species mutated to contain the reported SNPs passed quality assurance checks (Supplementary Table 4). This result may justify the use of only isotype 1 β-tubulins in in silico docking experiments. In addition, the alignment of both coding nucleotide and amino acid sequences revealed that in all species and isotypes, the positions holding the SNPs related to benzimidazole resistance had very similar codon constitutions, resulting in the same amino acids throughout (Fig. 1).

The cladograms derived from the phylogenetic analysis of nucleotide coding (Fig. 2A) and amino acid sequences (Fig. 2B) revealed that hookworm β-tubulins clustered based on their respective isotypes. Isotypes 1 and 2 in all species were more evolutionarily related, and the same was true for isotypes 3 and 4. As expected, the β-tubulins of N. americanus were genetically distant from those of the three Ancylostoma species based on the cladogram from the nucleotide sequences. Conversely, based on the tree derived from the amino acid sequences, A. duodenale β-tubulin isotype 4 was genetically distant from those of A. caninum, A. ceylanicum, and N. americanus. Moreover, based on the same tree, isotypes 3 and 4 showed greater evolutionary divergence than isotypes 1 and 2. The pairwise structural alignment of the predicted hookworm β-tubulin proteins revealed that all structure predictions had high structural similarity to the defined β-tubulin crystal structure from PDB ID No. 1jff (Supplementary Fig. 1).

3.2. SNPs affect benzimidazole binding affinities

In silico docking analysis revealed varying binding affinities between different hookworm species and benzimidazole drugs (Fig. 3). Among the β-tubulin isotype 1 proteins across all species, MBZ had the highest binding affinity, while ABZ had the lowest. Moreover, the binding affinities of the different benzimidazole ligands were lowest in N. americanus; the decrease was particularly marked for FBZ, OBZ, and ABZ reflecting the results of the phylogenetic analysis. Mutations of the β-tubulins through the introduction of the reported SNPs resulted in marked alterations in binding affinities. Overall, the Q134H and E198K mutations markedly decreased the binding affinities of the hookworm β-tubulins, regardless of the interacting benzimidazole drug. In addition, a similar but less reduction in binding affinity was observed when β-tubulins had the F200L mutation; the opposite was observed in proteins with the F200Y mutation. In contrast, the other mutations caused either slight increases or unchanged docking binding affinities.

3.3. SNPs affect benzimidazole binding positions and interacting residues

In addition to the altered binding affinities, SNPs also affected the binding pose between the β-tubulin proteins and benzimidazole drugs (Fig. 4). The observed reduction in binding affinities of hookworm β-tubulins with E198K and Q134H was supported by alterations in the docking poses of the protein and benzimidazole ligands. The loss of key carbon-hydrogen bonds was observed when glutamine was swapped with histidine at position 134. However, this mutation resulted in the same interaction as the wild-type counterpart. Taken together, these examples and the rest of the results show that resistance-associated mutations in the β-tubulin isotype 1 protein of different hookworm species caused alterations in binding positions when compared to ligand interactions with wild-type proteins. Moreover, changes in the types of bonds formed were also induced by the mutations. In particular, replacing glutamic acid with lysine at position 198 resulted in the emergence of Pi-cation bonds instead of hydrogen bonds. Lastly, shifts in binding positions due to mutations also affect interactions with other amino acids. For N. americanus β-tubulin with the E198K mutation, the mutant lysine formed conventional hydrogen bonds that reduced the drug interaction of phenylalanine at position 200 to noncovalent pi bonds. These changes in binding positions may result in the emergence of unfavorable interactions, changes in the type of bonds, and loss of key interactions thereby similarly affecting binding affinities.

The alteration of binding positions and interactions also resulted in several changes in the amino acids that interacted with benzimidazole drugs (Supplementary Tables 5 and 6). Leucine at position 253, methionine at 257, and alanine at position 314 were the most prevalent amino acid residues observed in all docking experiments with both wild-type and mutated proteins. Both wild-type (i.e., Glutamic Acid) and mutated (i.e., Alanine, Lysine, and Valine) amino acids at positions 198 and 134 were also observed to interact frequently with the drugs. These results revealed that among the positions that contain resistance-associated mutations, only those at positions 198 and 134 showed a significant number of interactions with ABZ, FBZ, MBZ, and OBZ collectively. These observations further strengthen the observed deviations in docking binding affinities, implicating mutations at positions 198 and 134 as markers of benzimidazole resistance.

3.4. E198K & Q134 mutations cause weaker benzimidazole binding dynamics

To further investigate the results of the molecular docking experiments, hookworm β-tubulin with the E198K and Q134H mutations complexed with FBZ and ABZ were analyzed using molecular dynamics simulations. Overall, the hookworm β-tubulin, both wild-type and mutated, and benzimidazole ligand complexes stabilized around 20 to 30 ns (Fig. 5A and B). The wildtype-, E198K- and Q134H-ABZ complexes have comparable RMSD values (wildtype = 2.32 ± 0.22, E198K = 2.30 ± 0.18, Q134H = 2.18 ± 0.26). Meanwhile, E198K and Q134H β-tubulin complexed with FBZ caused slight increases in RMSD when compared to the wildtype (wildtype = 1.99 ± 0.21, E198K = 2.37 ± 0.30, Q134H = 2.19 ± 0.17) over the 100 ns simulation. The RMSF values showed that the complexes shared a similar profile of flexible regions, regardless of the mutations. Among the β-tubulins complexed with ABZ, flexible regions included residues 37–45, 94–99, and 274–284 (Fig. 5C). For those docked with FBZ, increased flexibility was observed at residues 172–180, 244–253, and 275–282 (Fig. 5D). Both the RMSF and RMSD plots reveal that the wild-type and mutated β-tubulin-benzimidazole complexes share similar dynamic behaviors.

Hydrogen bonds between the β-tubulin receptors, both wild-type and mutated, and the benzimidazole drug ligands were determined during molecular dynamics simulations (Fig. 6). Overall, more H-bond contacts were formed between the wild-type and ABZ and FBZ. A decrease in the number of H-bonds was observed when the β-tubulin containing the E198K mutation interacted with ABZ. This decrease was more pronounced when the β-tubulin had the Q134H mutation. Using FBZ as a ligand, marked reductions in the number of H-bonds were observed when E198K and Q134H mutations were present in the β-tubulin receptor. These results indicate that the E198K and Q134H mutations cause weaker binding interactions with benzimidazole drugs indicating their role in benzimidazole resistance among hookworms.

The binding free energy is used to estimate the energetic differences between the bound and unbound states of the protein-ligand complexes, which can be used to quantify the stability and favorability of binding. In this study, the MM-PBSA and MM-GBSA methods were used to estimate the binding free energy between the β-tubulins and benzimidazole ligands (Supplementary Tables 7 and 8). The results showed that complexes with E198K and Q134H mutations in β-tubulins showed reduced binding energies compared to those with the wild-type (Fig. 7). Regardless of the benzimidazole ligand and method of binding free energy determination, the Q134H mutation caused significant decreases in binding energies, thereby causing weaker interactions. Taken together with results depicted in Fig. 6, the reduction in binding free energies among the complexes with the β-tubulin mutations is caused by diminished number of hydrogen bonds between the protein and the ligand. Molecular dynamics simulations proved that the E198K and Q134H mutations in hookworm β-tubulins confer benzimidazole resistance.

Benzimidazole resistance is an emerging threat that can potentially undermine efforts to control and eliminate STH, including hookworms. Despite concerns raised recently, the exact mechanism of resistance among these parasites remains unclear ^15,36. In this study, we used in silico docking with open-source platforms to study the potential consequences of SNPs in the β-tubulin isotype 1 gene of several hookworm species. Phylogenetic analysis revealed that the hookworm β-tubulin nucleotide and amino acid sequences clustered within their specific isotypes. In silico docking analysis revealed that amino acid substitutions brought about by SNPs in 4 codon locations of the hookworm β-tubulin isotype 1 protein resulted in altered binding affinities and binding positions. Among these mutations, E198K and Q134H may be important in conferring benzimidazole resistance. Molecular dynamics analysis revealed that the inclusion of these mutations did not significantly alter the binding dynamics of the β-tubulin-benzimidazole complex. However, the E198K and Q134H mutations caused marked reductions in the H-bonds formed, which resulted in significant reductions in the binding free energy. These results indicated that these mutations confer benzimidazole resistance.

The clustering of various β-tubulin genes based on their isotypes observed in this study was similar to those reported previously. Reconstruction of Ascaridomorpha phylogenies based on β-tubulin genes showed that ascarid β-tubulins were clearly separated into seven known isotypes (isotypes A–G) ²⁶. This clustering is believed to have arisen early in the evolution of various roundworm species, such as Ascaris lumbricoides, Anisakis simplex, and Parascaris equorum, within the Ascaridomorpha infraorder ²⁶. Our results, which indicated that hookworm β-tubulin isotypes 1 and 2 are monophyletic but distant from isotypes 3 and 4, which are also closely related, are also echoed by the results of Jones et al. ²⁵. They observed a similar clustering of strongyle β-tubulins and revealed a similar phylogenetic profile for other important parasitic nematodes, such as Haemonchus contortus and Trichuris trichiura. The results of our phylogenetic analysis and the aforementioned results from other published studies may have implications for benzimidazole resistance. Since benzimidazole resistance is often linked to SNPs in the β-tubulin isotype 1 gene of nematodes ^16,18, its close monophyletic relationship with isotype 2 suggests that SNPs occurring in the latter may also confer or assist in the emergence of benzimidazole resistance. This was previously observed by Doyle et al. ³⁹ when they studied resistant H. contortus. They observed that increasing Thiabendazole EC₅₀ concentrations correlated well with the E198V variant frequency in β-tubulin isotype 2 genes. Moreover, they noted that this mutation in isotype 2 may mediate a higher level of resistance than those conferred by SNPs in isotype 1 (e.g., F200Y). However, it is important to note that there are differences in expression levels of different nematode β-tubulins isotypes based on life stage and cell type ⁴⁰.

In silico docking of wild-type and mutated β-tubulin isotype 1 proteins with various benzimidazole drug ligands in this study revealed that resistance-associated mutations alter binding affinities and positions. Our results showed that E198K and Q134H mutations caused a marked decrease in binding affinity, regardless of the benzimidazole ligand used. The Q134H mutation was recently found to be widespread among A. caninum in the United States ¹⁹. The results of our study suggest that the same mutation may also confer resistance in A. ceylanicum, A. duodenale, and N. americanus. The E198K mutation, which caused a marked reduction in binding affinities in all hookworm species studied herein, has only been reported recently in A. caninum from Australia ²². The same SNP was recently reported in H. contortus infecting goats from several districts of Uganda ⁴¹. Fungal pathogens of crops in China that were resistant to benzimidazole derivatives also carry these mutations ⁴². Currently, there is no information regarding their effect on resistance among parasitic nematodes, which is a public health concern. Hence, for the first time, we provide evidence of marked binding affinity reduction caused by the E198K mutation in hookworm β-tubulin isotype 1 proteins with several benzimidazole drugs. Moreover, we also showed that the Q134H mutation may confer resistance in many hookworm species, as proven previously in A. caninum ¹⁹.

In the present study, we observed that mutations in the β-tubulin isotype 1 protein caused a shift in the ligand-binding position, modifications in bonds formed, and loss of important interactions. These alterations have been reported in previous in silico docking studies of other nematode species ^25,26. Despite the modifications in binding positions in mutated proteins, the profile of interacting amino acids with the benzimidazole drug ligands in this study remained similar to those of the wild-type proteins. Numerous drug interactions with leucine at position 253, methionine at position 257, and alanine at position 314 have also been noted in previous prediction and molecular dynamics studies ^35,43. Likewise, the frequent interactions with amino acids at 198, both in wild-type and mutated β-tubulins, observed in this study were also confirmed by prior research, thereby implicating that SNPs occurring at this position may be crucial in inducing resistance ^23,25,44. Our results, together with those previously reported, indicate that amino acid substitutions at positions 198 and 134 caused by resistance-related SNPs induce modifications in bond formation and other ligand-protein interactions that play roles in benzimidazole resistance in helminths of public and veterinary health significance. This should be explored further in future laboratory and field studies.

Molecular dynamics simulations revealed that hookworm β-tubulins, both wild-type and mutated, bound with benzimidazole ligands, share similar behavior in terms of atomic displacements and protein backbone flexibility. This indicated that resistance-associated mutations significantly altered the simulated behavior of the β-tubulin-benzimidazole ligand complex. This observation was also noted by Jones et al. ²⁶ in dynamics simulations of A. lumbricoides β-tubulin with the canonical mutations (i.e., F167Y, E198A, and F200Y) bound with ABZ. Despite these dynamic similarities, analysis of the H-bonds formed between the receptor and the benzimidazole ligand provided insights into how the E198K and Q134H mutations cause resistance. The two mutations decreased the number of H-bonds between the complexes. Benzimidazole resistance-associated mutations causing reduction in interactions have also been reported in in silico studies that assessed H. conturtos and Giargia lamblia β-tubulins, where residues cause intra-protein interactions instead of the ligand ⁴³. The marked differences in absolute binding free energy between complexes with wild-type and mutated β-tubulins also confirm the effects of lowering the number of H-bonds observed in this study. The reduced number of hydrogen bonds resulting in diminished binding free energy is due to the loss of attraction forces and structural stability, as observed in the analysis of different drug candidates docked with the hACE2-S protein ⁴⁵. Our molecular dynamics results provide an in silico background of how these mutations cause phenotypic resistance. Using C. elegans that were edited to have the Q134H mutation, which is widespread among A. caninum from the US, phenotypic resistance as measured by regressed animal length was proven after exposure to 30 µM of ABZ ¹⁹. Moreover, phenotypic resistance was also observed in C. elegans edited to contain E198K mutation when challenged with 30 µM of ABZ and FBZ ⁴⁶. Altogether, the molecular dynamics analysis of the β-tubulin-benzimidazole ligand complexes proved that the E198K and Q134H confer benzimidazole resistance by negatively affecting binding efficiency through reducing the number of interactions which drives down the absolute binding free energy.

Our results provided several insights into benzimidazole resistance in hookworms. Our phylogenetic analysis revealed that the hookworm β-tubulins clustered with respect to their isotypes. Hence, resistance-causing mutations in a particular β-tubulin isotype in one species may potentially cause the same as in other closely related species. Moreover, our in silico docking studies found that various amino acid substitutions brought about by resistance-related SNPs caused alterations in binding affinities and positions. In particular, E198K and Q134H caused marked reductions in binding affinities regardless of the hookworm species and interacting benzimidazole ligand. These mutations also resulted in shifts in the ligand-binding position, loss of interactions, and/or changes in the types of bonds formed. The molecular dynamics analysis of the β-tubulin-benzimidazole ligand complexes proved that the E198K and Q134H confer benzimidazole resistance by negatively affecting binding efficiency by reducing the number of interactions that drive down the absolute binding free energy. Our results support the notion that these SNPs occurring in the β-tubulin isotype 1 gene are implicated in causing resistance.

Supplementary material

The supplementary material for this article can be found at [DOI].

Ethical standards

This research does not require ethical approval for it did not use animal or human subjects.

Data availability

Data generated in this study can be made available upon request from the corresponding author.

World Health Organization (2023) Soil-transmitted helminth infections. https://www.who.int/news-room/fact-sheets/detail/soil-transmitted-helminth-infections. Published 2023. Accessed September 16
Traversa D (2012) Pet roundworms and hookworms: A continuing need for global worming. Parasit Vectors 5:91. 10.1186/1756-3305-5-91
Clements ACA, Alene KA (2022) Global distribution of human hookworm species and differences in their morbidity effects: a systematic review. Lancet Microbe 3(1):e72–e79. 10.1016/S2666-5247(21)00181-6
Colella V, Bradbury R, Traub R (2021) Ancylostoma ceylanicum. Trends Parasitol 37(9):844–845. 10.1016/j.pt.2021.04.013
Kyu HH, Abate D, Abate KH et al (2018) Global, regional, and national disability-adjusted life-years (DALYs) for 359 diseases and injuries and healthy life expectancy (HALE) for 195 countries and territories, 1990–2017: a systematic analysis for the Global Burden of Disease Study 2017. Lancet 392(10159):1859–1922. 10.1016/S0140-6736(18)32335-3
Jourdan PM, Lamberton PHL, Fenwick A, Addiss DG (2018) Soil-transmitted helminth infections. Lancet 391(10117):252–265. 10.1016/S0140-6736(17)31930-X
James SL, Abate D, Abate KH et al (2018) Global, regional, and national incidence, prevalence, and years lived with disability for 354 diseases and injuries for 195 countries and territories, 1990–2017: a systematic analysis for the Global Burden of Disease Study 2017. Lancet 392(10159):1789–1858. 10.1016/S0140-6736(18)32279-7
Stracke K, Jex AR, Traub RJ (2020) Zoonotic Ancylostomiasis: An Update of a Continually Neglected Zoonosis. Am J Trop Med Hyg 103(1):64–68. 10.4269/ajtmh.20-0060
Zibaei M, Nosrati MRC, Shadnoosh F et al (2020) Insights into hookworm prevalence in Asia: a systematic review and meta-analysis. Trans R Soc Trop Med Hyg 114(3):141–154. 10.1093/trstmh/trz115
Mationg MLS, Tallo VL, Williams GM et al (2021) The control of soil-transmitted helminthiases in the Philippines: the story continues. Infect Dis Poverty 10(1):85. 10.1186/s40249-021-00870-z
World Health Organization. Guideline: Preventive Chemotherapy to Control Soil-Transmitted Helminth Infections in at-Risk Population Groups. Geneva: World Health Organization (2017) https://apps.who.int/iris/handle/10665/258983. Accessed August 24, 2023
TroCCAP. Guidelines for the Diagnosis, Treatment and Control of Canine Endoparasites in the Tropics. 2nd ed (2019) https://www.troccap.com/2017press/wp-content/uploads/2019/05/TroCCAP_Canine_Endo_Guidelines_English_Ver2.pdf. Accessed August 23, 2023
TroCCAP. Guidelines for the Diagnosis, Treatment and Control of Feline Endoparasites in the Tropics. 2nd ed (2019) https://www.troccap.com/2017press/wp-content/uploads/2019/06/TroCCAP_Feline_Endo_Guidelines_English_Ver2.pdf. Accessed September 15, 2023
Montresor A, Mupfasoni D, Mikhailov A et al (2020) The global progress of soil-transmitted helminthiases control in 2020 and World Health Organization targets for 2030. Babu S. ed PLoS Negl Trop Dis 14(8):e0008505. 10.1371/journal.pntd.0008505
Tinkler SH (2020) Preventive chemotherapy and anthelmintic resistance of soil-transmitted helminths – Can we learn nothing from veterinary medicine? One Health 9:100106. 10.1016/j.onehlt.2019.100106
Von Samson-Himmelstjerna G, Blackhall WJ, McCARTHY JS, Skuce PJ (2007) Single nucleotide polymorphism (SNP) markers for benzimidazole resistance in veterinary nematodes. Parasitology 134(8):1077–1086. 10.1017/S0031182007000054
Prichard RK (2007) Markers for benzimidazole resistance in human parasitic nematodes? Parasitology 134(8):1087–1092. 10.1017/S003118200700008X
Furtado LFV, De Paiva Bello ACP, Rabelo ÉML (2016) Benzimidazole resistance in helminths: From problem to diagnosis. Acta Trop 162:95–102. 10.1016/j.actatropica.2016.06.021
Venkatesan A, Castro PDJ, Morosetti A et al (2023) Molecular evidence of widespread benzimidazole drug resistance in Ancylostoma caninum from domestic dogs throughout the USA and discovery of a novel β-tubulin benzimidazole resistance mutation. PLoS Pathog 19(3):e1011146. 10.1371/journal.ppat.1011146
Evason MD, Weese JS, Polansky B, Leutenegger CM (2023) Emergence of canine hookworm treatment resistance: Novel detection of Ancylostoma caninum anthelmintic resistance markers by fecal PCR in 11 dogs from Canada. Am J Vet Res 84(9). 10.2460/ajvr.23.05.0116
Castro PDJ, Durrence K, Durrence S et al (2023) Multiple anthelmintic drug resistance in hookworms (Ancylostoma caninum) in a Labrador breeding and training kennel in Georgia, USA. J Am Vet Med Assoc 261(3):342–347. 10.2460/javma.22.08.0377
Stocker T, Scott I, Šlapeta J Unambiguous identification of Ancylostoma caninum and Uncinaria stenocephala in Australian and New Zealand dogs from faecal samples. Aust Veterinary J. August 2023:avj.13272. 10.1111/avj.13272
Medeiros CDS, Furtado LFV, Miranda GS, Da Silva VJ, Dos Santos TR, Rabelo ÉML (2022) Moving beyond the state of the art of understanding resistance mechanisms in hookworms: confirming old and suggesting new associated SNPs. Acta Trop 233:106533. 10.1016/j.actatropica.2022.106533
Furtado LFV, de Aguiar PHN, Zuccherato LW et al (2019) Albendazole resistance induced in Ancylostoma ceylanicum is not due to single-nucleotide polymorphisms (SNPs) at codons 167, 198, or 200 of the beta-tubulin gene, indicating another resistance mechanism. Parasitol Res 118(3):837–849. 10.1007/s00436-019-06218-9
Jones BP, van Vliet AHM, LaCourse EJ, Betson M (2022) In Silico Docking of Nematode β-Tubulins With Benzimidazoles Points to Gene Expression and Orthologue Variation as Factors in Anthelmintic Resistance. Frontiers in Tropical Diseases. ;3. https://www.frontiersin.org/articles/10.3389/fitd.2022.898814. Accessed October 26, 2023
Jones BP, van Vliet AHM, LaCourse EJ, Betson M (2022) Identification of key interactions of benzimidazole resistance-associated amino acid mutations in Ascaris β-tubulins by molecular docking simulations. Sci Rep 12(1):13725. 10.1038/s41598-022-16765-4
Dantas-Torres F, Ketzis J, Mihalca AD et al (2020) TroCCAP recommendations for the diagnosis, prevention and treatment of parasitic infections in dogs and cats in the tropics. Vet Parasitol 283:109167. 10.1016/j.vetpar.2020.109167
Howe KL, Bolt BJ, Shafie M, Kersey P, Berriman M (2017) WormBase ParaSite – a comprehensive resource for helminth genomics. Mol Biochem Parasitol 215:2–10. 10.1016/j.molbiopara.2016.11.005
The UniProt Consortium (2019) UniProt: a worldwide hub of protein knowledge. Nucleic Acids Res 47(D1):D506–D515. 10.1093/nar/gky1049
Hall TA (1999) BioEdit: a user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT. Nucleic acids symposium series. ;41(41):95–98
Thompson JD, Higgins DG, Gibson TJ (1994) CLUSTAL W: improving the sensitivity of progressive multiple sequence alignment through sequence weighting, position-specific gap penalties and weight matrix choice. Nucleic Acids Res 22(22):4673–4680. 10.1093/nar/22.22.4673
Tamura K, Stecher G, Kumar S (2021) MEGA11: Molecular Evolutionary Genetics Analysis Version 11. Mol Biol Evol 38(7):3022–3027. 10.1093/molbev/msab120
Schwede T, Kopp J, Guex N, Peitsch MC (2003) SWISS-MODEL: an automated protein homology-modeling server. Nucleic Acids Res 31(13):3381–3385. 10.1093/nar/gkg520
Bueno O, Estévez Gallego J, Martins S et al (2018) High-affinity ligands of the colchicine domain in tubulin based on a structure-guided design. Sci Rep 8(1):4242. 10.1038/s41598-018-22382-x
Robinson MW, McFerran N, Trudgett A, Hoey L, Fairweather I (2004) A possible model of benzimidazole binding to β-tubulin disclosed by invoking an inter-domain movement. J Mol Graph Model 23(3):275–284. 10.1016/j.jmgm.2004.08.001
Ng’etich AI, Amoah ID, Bux F, Kumari S (2023) Anthelmintic resistance in soil-transmitted helminths: One-Health considerations. Parasitol Res 123(1):62. 10.1007/s00436-023-08088-8
Trott O, Olson AJ, AutoDock, Vina (2010) Improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading. J Comput Chem 31(2):455–461. 10.1002/jcc.21334
Kim S, Thiessen PA, Bolton EE et al (2016) PubChem Substance and Compound databases. Nucleic Acids Res 44(D1):D1202–D1213. 10.1093/nar/gkv951
Doyle SR, Laing R, Bartley D et al (2022) Genomic landscape of drug response reveals mediators of anthelmintic resistance. Cell Rep 41(3). 10.1016/j.celrep.2022.111522
Saunders GI, Wasmuth JD, Beech R et al (2013) Characterization and comparative analysis of the complete Haemonchus contortus β-tubulin gene family and implications for benzimidazole resistance in strongylid nematodes. Int J Parasitol 43(6):465–475. 10.1016/j.ijpara.2012.12.011
Kalule F, Vudriko P, Nanteza A et al (2023) Prevalence of gastrointestinal parasites and molecular identification of beta-tubulin mutations associated with benzimidazole resistance in Haemonchus contortus in goats from selected districts of Uganda. Veterinary Parasitology: Reg Stud Rep 42:100889. 10.1016/j.vprsr.2023.100889
Liu Y, Chen X, Jiang J, Hamada MS, Yin Y, Ma Z (2014) Detection and dynamics of different carbendazim-resistance conferring β-tubulin variants of Gibberella zeae collected from infected wheat heads and rice stubble in China. Pest Manag Sci 70(8):1228–1236. 10.1002/ps.3680
Aguayo-Ortiz R, Méndez-Lucio O, Romo-Mancillas A et al (2013) Molecular basis for benzimidazole resistance from a novel β-tubulin binding site model. J Mol Graph Model 45:26–37. 10.1016/j.jmgm.2013.07.008
Aguayo-Ortiz R, Méndez-Lucio O, Medina-Franco JL et al (2013) Towards the identification of the binding site of benzimidazoles to β-tubulin of Trichinella spiralis: Insights from computational and experimental data. J Mol Graph Model 41:12–19. 10.1016/j.jmgm.2013.01.007
Pandey P, Rane JS, Chatterjee A et al (2021) Targeting SARS-CoV-2 spike protein of COVID-19 with naturally occurring phytochemicals: an in silico study for drug development. J Biomol Struct Dynamics 39(16):6306–6316. 10.1080/07391102.2020.1796811
Dilks CM, Koury EJ, Buchanan CM, Andersen EC (2021) Newly identified parasitic nematode beta-tubulin alleles confer resistance to benzimidazoles. Int J Parasitology: Drugs Drug Resist 17:168–175. 10.1016/j.ijpddr.2021.09.006

The authors declare no competing interests.

SupplementaryMaterials.docx

Download PDF

Version 1

posted

You are reading this latest preprint version

Unraveling the Mechanisms of Benzimidazole Resistance in Hookworms: A Molecular Docking and Dynamics Study

Status:

Version 1

Abstract

Background

Methods

Results

Conclusion

Figures

1. Introduction

2. Materials and Methods

2.1. Hookworm β-tubulin nucleotide and amino acid sequences

2.2. Phylogenetic Analysis

2.3. Homology Modelling and Quality Assurance Checks

2.4. In silico Docking

2.5. Molecular dynamics simulations

3. Results

3.1. Hookworm β-tubulins and their phylogeny

3.2. SNPs affect benzimidazole binding affinities

3.3. SNPs affect benzimidazole binding positions and interacting residues

3.4. E198K & Q134 mutations cause weaker benzimidazole binding dynamics

4. Discussion

5. Conclusion

Declarations

Supplementary material

Ethical standards

Data availability

References

Additional Declarations

Supplementary Files

Status:

Version 1