In Silico Docking of Nematode β-Tubulins With Benzimidazoles Points to Gene Expression and Orthologue Variation as Factors in Anthelmintic Resistance

The efficacy of benzimidazole anthelmintics can vary depending on the target parasite, with Ascaris nematodes being highly responsive, and whipworms being less responsive. Anthelmintic resistance has become widespread, particularly in strongyle nematodes such as Haemonchus contortus in ruminants, and resistance has recently been detected in hookworms of humans and dogs. Past work has shown that there are multiple β-tubulin isotypes in helminths, yet only a few of these contribute to benzimidazole interactions and resistance. The β-tubulin isotypes of ascarids and soil-transmitted helminths were identified by mining available genome data, and phylogenetic analysis showed that the ascarids share a similar repertoire of seven β-tubulin isotypes. Strongyles also have a consistent pattern of four β-tubulin isotypes. In contrast, the whipworms only have two isotypes, with one of these clustering more basally and distinct from any other group. Key β-tubulin isotypes selected based on previous studies were the focus of in silico molecular docking simulations to look at the interactions with benzimidazoles. These showed that all β-tubulins had similar interactions with benzimidazoles and maintained the key bond with residue E198 in all species, indicating similar mechanisms of action. However, the interaction was stronger and more consistent in the strongyles and whipworms than it was in the ascarids. Alteration of β-tubulin isotypes with the common resistance-associated mutations originally identified in H. contortus resulted in similar interaction modeling for all species. In conclusion, ascarids, strongyles, and whipworms all have their own unique repertoire of β-tubulins, which could explain why benzimidazole resistance and susceptibility varies between these groups of parasites. These data complement recent work that has highlighted the roles of essential residues in benzimidazole drug binding and shows that there is a separation between strongyle parasites that frequently develop resistance and ascarid parasites, which have been much less prone to developing resistance.


INTRODUCTION
Due to the early emergence of benzimidazole (BZ) resistance in nematodes of veterinary importance, such as Haemonchus contortus, and its link with b-tubulins, there is growing interest in the diversity of b-tubulins in nematodes (1,2). The model organism Caenorhabditis elegans, a free-living nematode, is the best source of information regarding the diversity of b-tubulins and the mechanisms of resistance. There are six b-tubulin genes (isotypes) within the genome of C. elegans and other helminths (1,(3)(4)(5)(6)(7)(8)(9). Benzimidazole resistance was specifically associated with one of the isotypes called ben-1 in C. elegans, and orthologues of this isotype have been associated with BZ resistance in strongyle helminths of veterinary importance, such as Cooperia oncophora, Haemonchus contortus, and Ostertagia ostertagia (1,4). There are currently four known b-tubulin isotypes in H. contortus, and non-strongyle parasites such as Trichuris and Ascaris also have multiple b-tubulins, although in these helminths, the b-tubulin isotypes linked with BZ interactions are not orthologues of C. elegans ben-1 (1,4).
Treatment of the ruminant nematode H. contortus has been most impacted by drug resistance and has subsequently been the focus of most research. Previous work has shown that of the four b-tubulins found in the H. contortus genome, only isotypes 1 and 2 are expressed in adults at high enough levels to be impacted by BZs, and isotype 1 has been the focus of research into drug resistance (1). A range of mutations in three amino acids have been shown to be associated with BZ resistance in several nematode species (10)(11)(12)(13)(14)(15)(16)(17). A limited number of studies have used in silico methods to identify the exact mechanisms at play in b-tubulin-BZ interactions and the effects of resistance-associated mutations on these interactions (18,19).
Recent works in Ascaris have provided the first in-depth investigations into the b-tubulins of Ascaris, their role in BZ interactions, and the potential consequences of resistanceassociated mutations (7,20). We now know that similar to strongyle parasites, Ascaris has several highly expressed isotypes capable of interacting with BZs, with isotype A being the most likely candidate (7). In silico work has highlighted the differences in drug binding between wild-type and mutated Ascaris b-tubulins, providing a basis for future work on parasite b-tubulins (20). There have also been multiple studies on Parascaris and Ascaridia, both members of the ascarid family of helminths, which have shown that resistance is emerging in these ascarids, though the mutations responsible are yet to be confirmed (9,21,22).
The efficacy of BZs against the soil-transmitted helminth (STH) Ascaris has remained generally high. However, this is not the case for other STHs such as hookworms and whipworms. Meta-analysis found that while average cure rates for Ascaris lumbricoides were between 88% and 95%, for hookworm, they were lower: between 15% and 72% depending on which BZ was used. Whipworm showed consistently lower susceptibility to BZ with average cure between 28% and 36% (23). More recent studies have continued to find lower susceptibility of whipworms to BZ, suggesting that whipworms may be innately less susceptible (24)(25)(26). The typical resistance-associated single-nucleotide polymorphisms (SNPs) (F167Y, E198A, and F200Y) have been found in both hookworms and whipworms, demonstrating that there has likely been a selective pressure by BZs leading to reduced treatment response (27)(28)(29)(30)(31).
While there have been investigations into BZ docking in a number of organisms, including parasitic nematodes and BZsusceptible and refractory protozoan parasites and mammals, the recent work on Ascaris was the first to look into STHs of public health importance (18)(19)(20)32). No work has been carried out on the remaining two groups of STH, the hookworms (Ancylostoma and Necator americanus) and whipworms (Trichiura), even though BZ resistance is a known problem in the control of STHs.
The aim of this work was to identify and compare the relationships between the b-tubulins of STHs and Ascaridomorpha parasites, herein referred to as ascarids, and use those isotypes most likely to interact with BZs to perform in silico drug docking simulations. We aimed to determine if all ascarids share the same grouping of b-tubulins as Ascaris and see how this compares to other groups of helminths. We also aimed to investigate whether there are any intrinsic differences in the binding capacity of ascarid, hookworm, and whipworm btubulins that might explain the varying efficacies of BZs against different nematodes.

b-Tubulin Identification
The genomes of common STHs and all ascarid species with genomes accessible in Wormbase-Parasite were searched to identify b-tubulin isotypes (33,34). A BLAST search of the genomes was performed using the A. lumbricoides isotype A peptide sequence as a reference (35). Genomes were also searched using the term "Tubulin beta" to ensure as many genes as possible could be found. The results of this were then filtered based on size and tubulin type to remove small fragment sequences and atubulins. The published b-tubulins of ascarids available from NCBI were also included (36). The nematode species included in this analysis were as follows: Ancylostoma caninum, Ancylostoma ceylanicum, Ancylostoma duodenale, Necator americanus, Trichuris trichiura, Trichuris suis, Parascaris equorum, Parascaris univalens, Anasakis simplex, Toxocara canis, and Ascaridia galli. Although Asc. galli does not have a genome available, the sequence for a btubulin was found in the BLAST search. A full list of b-tubulin genes identified can be seen in Supplementary Table 1

Phylogenetic Analysis
To compare the genetic structuring of the identified b-tubulins, phylogenetic analysis of the peptide sequences was performed in MEGA X (38). The four b-tubulins described in H. contortus and the six b-tubulins from C. elegans were included in the phylogenies as a comparison of a better studied organism, and the b-tubulin sequences described previously in Ascaris were included to show the full b-tubulin repertoire of ascarids (7,20). The sequences were aligned using the MUSCLE algorithm, the best evolutionary model was found, and a maximum likelihood tree was constructed using the LG+G model with 1,000 bootstraps.

In Silico Simulations
Molecular dynamics simulations were carried out on the b-tubulin isotype 1 proteins for the hookworms and whipworms, and isotype A proteins for the ascarids. Isotype A was chosen for this analysis as it is the isotype that is associated with the highest expression levels in ascarids and is therefore the most likely to be affected by BZs. For the other STHs, isotype 1 was chosen as this is the one that has previously been shown to cause resistance when mutated. To create 3D structures of the b-tubulin sequences, homology models were constructed in SWISS-MODEL using the PDB structure 6fkj as a reference (available at: https://swissmodel. expasy.org/ [Accessed February 22, 2021]) (39)(40)(41). All homology models were minimized using the YASARA energy minimization server (available at: http://www.yasara.org/minimizationserver. htm [Accessed October 6, 2021]) and quality checks of the structures were carried out using ProSA-web (available at: https://prosa.services.came.sbg.ac.at/prosa.php [Accessed October 6, 2021]), Verify3D, and PROCHECK (42)(43)(44)(45)(46)(47). Verify3D and PROCHECK are part of the UCLA SAVES v6.0 server (available at: https://saves.mbi.ucla.edu/ [Accessed October 6, 2021]). Default settings were used for all servers.
Albendazole sulfoxide (ABZSO) was used to perform molecular dynamics simulation using Molecular Operating Environment (MOE) version 2020.01 (48). This drug was chosen as it is the active metabolite of albendazole, one of the most commonly administered BZs to treat helminths. Molecular dynamics simulations were carried out in MOE following the methods set out in Jones et al. (20). Briefly, b-tubulin structures were optimized using the Protonate3D command within MOE, the drug binding pockets were identified using the Site Finder algorithm, and ABZSO was docked into the pocket corresponding to the known binding region of BZs. The 10 best docking poses were returned based on the GBVI/WSA dG scoring method, and from this, the best binding pose was chosen for molecular dynamics simulations. Molecular dynamics simulations were carried out using the NPA algorithm and the Amber10:EHT forcefield. An equilibration step of 100 picoseconds (ps) at 300°K was run initially, followed by a further 500-ps run with a time step of 0.002 ps. Samples were taken for comparison every 50 ps and these structures were minimized to find drug binding affinity; the pose with the strongest affinity was taken as the final result.
The three canonical resistance-associated mutations (F167Y, E198A, and F200Y) were introduced into the Anc. duodenale and Tr. trichiura sequences by manually changing the sequence residues. Homology models were created and minimized and molecular dynamics simulations were carried out using the same methods as previously stated.

Nematode b-Tubulins Cluster Into Different Families
The initial genome search resulted in the identification of 107 btubulin sequences in 13 genomes, representing An. simplex, P. equorum, P. univalens, T. canis, Anc. caninum, Anc. ceylanicum, Anc. duodenale, N. americanus, Tr. trichiura, and Tr. suis. Additional sequences for Asc. galli and P. equorum were added from the NCBI database. Deduplication and removal of most incomplete sequences and a-tubulins left 64 b-tubulin sequences, which included the 14 Ascaris and four H. contortus reference sequences (Supplementary Table 1). Some incomplete sequences were included in the analysis as long as they were more than 200 amino acids long and included the region of the protein that contained the predicted binding pocket and resistance-associated SNP sites in the corresponding gene.
The phylogenetic analysis showed that the ascarid b-tubulins resembled the Ascaris b-tubulins, with representatives being present in each isotype cluster ( Figure 1). Parascaris univalens was the only organism containing all seven b-tubulin isotypes. Toxocara canis lacked isotype C or D, and for An. simplex, which had one of the lowest qualities of genome assembly according to the BUSCO and CEGMA metrics, isotypes D, E, F, and G could not be identified. This may well be a result of the lower quality of the available genomes, rather than a lack of these isotypes being present in these species. There were several additional isotypes found only in T. canis that were very similar to isotype A. Of these, two had large insertions and so were disregarded, but one was of a similar size to all other b-tubulin isotypes and so was designated as isotype A.
The non-ascarids (C. elegans, hookworms, H. contortus, and Trichuris) are more distantly related to Ascaris than the ascarids, and their b-tubulin repertoire reflected this. The strongyle genomes contained sequences that clustered with isotypes E and F, isotypes 4 and 3, respectively. However, within these isotype groupings, they formed their own clusters, with the corresponding C. elegans btubulins separate from the other ascarids. In the strongyles, isotypes 1 and 2 formed a distinct cluster separate from the ascarid isotypes. In C. elegans, isotypes 1 and 2 did not correspond to the strongyle btubulins, instead forming their own cluster closely related to isotype 1. The C. elegans ben-1 isotype was isolated basal to the strongyle isotypes 1 and 2 cluster. For Anc. duodenale, only isotype 1 was found, which may be due to the incompleteness of the genome ( Figure 1). There were some nodes that showed relatively low bootstrap values.
The Trichuris genomes contained only two b-tubulin isotypes. Isotype 1 clustered basal to the b-tubulin E/4 and F/3 clade. The second isotype did not cluster with any other group and was shown to be the most basal of the b-tubulins analyzed ( Figure 1). A phylogeny was also created omitting the divergent Trichuris isotype 2 sequences, which resulted in no major changes and so was discarded (data not shown).

Similarities in the Drug Interactions of Nematode b-Tubulins
Quality checks of the structures allow us to assess how realistic the homology models are compared to experimentally determined proteins by employing various tests. The quality checks of minimized homology models showed acceptable results with Zscores within the expected range and verify3D scores over 80% (Supplementary Figures 2-17, Supplementary Tables 2 and 3). Most models showed that all residues within the proteins were in allowable conformations as PROCHECK found no residues with unfavorable torsion angles; however, the Trichuris isotype 1 proteins showed some outlier residues in the PROCHECK Ramachandran plots (Supplementary Tables 2 and 3).
Molecular dynamics results showed that there were many interactions between the b-tubulins and BZs shared across species. All models showed strong binding between ABZSO and residue E198, except for P. equorum, where only a relatively weak bond between BZ and E198 was seen (Figures 2 and 3, Table 1). Interactions with Q134, L253, and N256 were the most frequently seen after E198, with bonds seen in four, five, and four b-tubulins, respectively. Simulations for isotypes from the model organism Caenorhabditis elegans, common soil-transmitted helminths, and other Ascaridomorpha species. The phylogeny was created using an LG+G model with 1,000 bootstraps. Ascaridomorpha isotypes were named based on Ascaris nomenclature and the strongyle and Trichuris isotypes were numbered based on established nomenclature. The Ascaridomorpha species follow a similar pattern of b-tubulin isotypes to Ascaris, organizing into seven groups [clusters (A-G) colored separately based on Ascaris isotype nomenclature]. In some cases, isotypes are not present for all species of Ascaridomorpha, likely due to the incompleteness of the available genomes. The strongyles (H. contortus and the hookworms) have their own b-tubulin profiles; however, these are closely related to the b-tubulins seen in the C. elegans. The strongyles contain four isotypes, two cluster distinctly from any Ascaridomorpha group (isotype 1 and 2), one clusters with group E (isotype 4), and one clusters with group F (isotype 3). The Trichuris have the most distinct b-tubulin profile, with only two isotypes found. One of these clusters basally to groups (E, F), and the other is distinct from any other b-tubulin isotype and forms its own separate group. Ascaris a-tubulins were used as the outgroup. the ascarid isotype A structures showed similar results to each other and to the molecular dynamic simulations performed for Ascaris described previously (20) (Figure 2, Table 1). The results for the hookworm and whipworm species showed that the strong bond with E198 was seen in all species (Figure 3, Table 1

Mutated Models Highlight Key Changes That May Affect Drug Binding
Molecular dynamics with the Anc. duodenale F167Y model showed that the strong bond between BZ and E198 was still present. A self-binding interaction was observed between Q134 and F167Y, and a strong bond with K350 was also seen ( Figure 4A, Table 2). In the Anc. duodenale E198A model, the bond between BZ and residue 198 was lost, and no bonds of similar strength were present ( Figure 4B, Table 2). The Anc. duodenale F200Y model showed that the E198 bond was not present, and a self-binding interaction was seen between E198 and F200Y. The strong bond with K350 was again present, The table shows the type of bonds formed (H-hydrogen bond, A-arene bond), the amino acid the bond is formed with, and the drug used (ABZSO = albendazole sulfoxide). The energy of the bonds between the amino acid and drug is given in kcal/mol, the distance between the bonded atoms is given in angstroms (Å), and the number of bonds formed between the amino acid and the drug is shown (frequency). The binding affinities of the protein-drug interactions from molecular dynamics simulations are shown in kcal/mol. The table shows that residue E198 (Glu198) is present in all models. The interaction energies with E198 can be seen to be noticeably higher that any of the other interactions, with the exception of K252 and K350, which only occasionally show high binding energy. The overall affinity of the protein-BZ binding is also shown. Due to a deletion at residue 177 in Trichuris isotype 1, residues 197, 199, and 255 correlate to 198, 200, and 256 in the other nematodes.
indicating that this residue may have an important role in drug stability in the binding pocket ( Figure 4C, Table 2). The results of molecular dynamics simulations for the Tr. trichiura F167Y model showed the presence of BZ binding with E197. The mutated F167Y residue also formed a bond with BZ ( Figure 4D, Table 2). In the E197A model, the key interaction with residue 197 was not present ( Figure 4E, Table 2). The residues that BZ did bind with meant that the drug was pulled away from the usual binding region and no interactions were made with residues 167, 197, and 199 (167, 198, and 200 in all other species). In the Tr. trichiura F199Y model, BZ binds with the key E197 residue. However, in this model, an ionic bond is formed rather than a hydrogen bond, which increases the binding energy ( Figure 4F, Table 2). A self-binding interaction was also seen between E197 and F199Y in Tr. trichiura ( Figure 4F).

DISCUSSION
This research has identified multiple b-tubulin isotypes in several parasites of human and veterinary importance. We have shown using phylogenetic analyses that the b-tubulin profiles are similar for the ascarids with seven groupings of b-tubulin isotypes. The strongyles also share similar b-tubulin profiles to each other using only four isotypes. The in silico work performed has shown that BZs appear to bind in the same way in all nematode b-tubulins investigated here and highlights the importance of residue 198 in the BZ mechanism of action. Due to the understudied nature of the organisms investigated here, in-depth genome assemblies are not available for all species. There were major disparities in the quality of the genomes used for this study, based on the quality metrics in Wormbase-Parasite and the number of tubulin genes that were recovered. For Asc. galli, a genome was not available, although a single sequence for the b-tubulin gene that had been previously studied in relation to BZ resistance was available (6). As Asc. galli is an important parasite of poultry, and with the recent report of BZ resistance being seen in a related species, Ascaridia dissimilis, it was determined that it would be important and of wider interest to include this single sequence in the study (22).
The clearest differences were seen in the Trichuris b-tubulins. Only two isotypes were identified: one that clusters basal to groups E and F of the other nematode b-tubulins, and another that clusters basally to the other helminth b-tubulins. This would suggest that Trichuris diverged from the other helminths prior to the b-tubulin duplications that gave rise to the b-tubulin families seen in the rest of the helminths. This fits with the taxonomy that shows that Trichuris are of a different class (Enoplea) to the strongyles and ascarids (Chromadorea). However, with this in mind, we should then expect the Trichuris isotype 1 to be basal to all Chromadorea b-tubulins, something that will require further investigation. In this study, we have used the current nomenclature for Trichuris b-tubulins of isotype 1 and 2. In the past, tubulins have been named based on the order in which they were discovered, or their expression level, and this has caused problems in the nomenclature, with the b-tubulin isotype 1 from different groups not being orthologues (1,4,8,9). This problem has been addressed recently for Ascaris where the naming convention was changed from a numerical system to alphabetical (A-G) to clarify the difference in evolutionary history and roles of b-tubulins in these groups of parasites (7). The Trichuris b-tubulin isotypes 1 and 2 are also not orthologues of the strongyle btubulins, and the argument that the naming of these isotypes is misleading could be made. However, the reclassification of the Trichuris b-tubulin nomenclature is beyond the scope of this study and would require further examination of other species more closely related to Trichuris to put forward a new nomenclature that would best fit these species. The full repertoire of b-tubulins characterized in Ascaris was also seen in P. univalens but none of the other ascarids. It is likely that this is the result of incomplete genome assemblies. This is evidenced by looking at the isotypes found in An. simplex and T. canis, which together represent all isotypes except isotype D. As Asc. galli is the only non-Ascaridoidea included in this group and a genome is not available, the only evidence to suggest that it would contain most, if not all the Ascaris isotypes is the fact that the strongyles, which are more distinct, also share some of these isotypes. Within the ascarids group, T. canis does show an extra isotype within group A. Without further investigation, it cannot be concluded if this is the result of recent duplication in this species or problems with genome assembly.
The phylogenies show that there are noticeable differences between the structuring of the ascarid and strongyle b-tubulins included in this study. There are some noticeable differences seen between the phylogeny presented here and phylogenies presented in recent studies of ascarid b-tubulins (7,9,49). While the general ordering of isotypes in our study is similar to that reported in previous studies, the exact clustering and evolutionary relationship is different. In our phylogeny, we see a continuous branching off of isotypes A, C, and D, with B and G, and E and F forming separate subsequent clusters. In previous work A, C, and D cluster together with B and G forming a separate cluster, and they together form a larger clade separate from E and F (7,9,49). Again, in previous work, the strongyle The table shows the type of bonds formed (H-hydrogen bond, IH-ionic hydrogen bond, A-arene bond), the amino acid the bond is formed with, and the drug used (ABZSO = albendazole sulfoxide). The energy of the bonds between the amino acid and drug is given in kcal/mol, the distance between the bonded atoms is given in angstroms (Å), and the number of bonds formed between the amino acid and the drug is shown (frequency). The binding affinities of the protein-drug interactions from molecular dynamics simulations are shown in kcal/ mol. Due to a deletion at residue 177 in Trichuris isotype 1, residues 197, 199, 235, 238, 246, 252, and 255 correlate to 198, 200, 236, 239, 247, 253, and 256 in the other nematodes. The tables show that both F167Y models retain the strong bond with E198 seen in the wild-type models; however, in Tr. trichiura, the mutated residue bonds with ABZSO. Both species E198A models show that this mutation leads to a loss of interaction between ABZSO and residue 198. In the F200Y models, ABZSO interactions with residue 198 return for Tr. trichiura but not for Anc. duodenale. The bond between E198 and ABZSO in the Tr. trichiura F200Y model has a much higher binding energy than wild-type models due to the ionic hydrogen bonding.
isotype 1 and 2 clusters were found to be at a basal position within the larger A, C, D, B, and G clade, whereas in our phylogeny, we see them form a separate clade (7,9,49). Within this clade, the relationship between the strongyles and C. elegans corresponds to previous work with C. elegans isotypes 1 and 2 being closely related to but isolated from strongyle isotype 1 and C. elegans ben-1 being basal to this clade. Trichuris isotype 1 is basal to isotypes E and F as was seen previously. However, in this former analysis, Tr. trichiura was also basal to all other nematode b-tubulins (7). Conversely, another investigation of nematode tubulins showed Tr. trichiura isotype 1 to sit within a clade with A. suum isotype E and C. elegans isotype 4 (4). Finally, the relationship between ascarid isotype F and strongyle isotype 3 and between ascarid isotype E and strongyle isotype 4 matches that seen previously with these isotypes forming their own clade (7,9,49). There are several possible causes for the differences seen between our phylogeny and those previously published such as the method used to construct the phylogeny, or the sequences used. Our study included more sequences than previous analyses have, including more divergent isotypes such as Trichuris isotype 2, C. elegans isotype 6, and Ascaris a-tubulin, which give us an insight into the history of b-tubulin evolution. The inclusion of these sequences may be the reason for the slightly different clustering of isotypes and the low bootstrap values seen for some nodes, although no major differences were seen when Trichuris isotype 2 was removed from analysis. The isotypes found in the hookworm genomes match the four b-tubulin isotypes that have been previously identified in H. contortus. Previous work has shown that of the four isotypes, only isotype 1 and isotype 2 are highly expressed in H. contortus (1). As most BZ resistance is linked to isotype 1 in strongyles, it was this isotype that was chosen to be the focus of further experiments.
In the molecular dynamics simulations performed in this study, most of the nematode b-tubulins show similar BZ binding results to the H. contortus and Ascaris models in previous studies (19,20). A strong bond formation is seen between BZ and the E198 residue, which is key to BZ sensitivity. The exceptions to this were in P. equorum, which had a weaker bond with E198 compared to other helminths. While, so far, there has been a lack of evidence for BZ resistance in Ascaris, in Parascaris and Ascaridia, reports of resistance are beginning to be seen (9,22,50,51). However, there is no molecular evidence so far of the canonical BZ resistance-associated SNPs in resistant populations of ascarids (9).
The hookworms have the same repertoire of b-tubulins as the more well-studied strongyles such as H. contortus. Most importantly, they share the two closely related isotypes associated with BZ interaction that are highly expressed (1). The theory of redundancy has been posited for b-tubulins in C. elegans, whereby the loss of function of either b-tubulin 1 or 2 rarely causes cellular defects whereas missense mutations lead to more severe effects (3,52). Recent work has also proven that the canonical resistance-associated mutations cause resistance in gene-edited C. elegans b-tubulins yet have no effect on fitness in controls (53). It could be the case that as there are two proteins performing the same job in strongyles, this allows one to gain mutations that confer resistance but may reduce fitness, while the other protein can take on the needs of the cells counteracting this burden.
Most evidence of BZ resistance has come from strongyle parasites with mutations in isotypes 1 and 2 being the main cause (1,10,12,13,(54)(55)(56)(57). These two isotypes do not group closely with any ascarid b-tubulin isotype (Figure 1). It may be the close genetic relationship between these isotypes that provides the redundancy needed for mutations to develop. The ascarids have the most b-tubulin isotypes, yet BZs still seem to work efficiently against them. The greater divergences between the ascarid btubulin isotypes may make it less likely that they can take on the roles of isotype A and therefore the fitness cost of the resistanceassociated SNPs would make them deleterious.
Trichuris are much more genetically divergent from both ascarids and strongyles, and hence have a very distinct repertoire of b-tubulins. Although there is little difference seen in terms of overall binding and affinity between Trichuris b-tubulins and BZs, the mechanisms may be different in these parasites. For example, it is possible that both b-tubulins are highly expressed, explaining why normal doses of BZs have been ineffective, yet higher doses or prolonged treatment increases the efficacy of BZs in Trichuris (25,58).
Although there has been a documented decline in the efficacy of BZs against whipworms, there is still little evidence for the presence of the canonical resistance-associated SNPs (59). In several studies, these SNPs could not be found; in one study, only the F200Y SNP was found, and in another, all three SNPs were found in different populations (27,(60)(61)(62)(63)(64). In a population where all three SNPs were found, the proportion of the F167Y SNP was decreased after treatment (27). While there has been a slower development of resistance in hookworms compared to ruminant strongyles, reports are beginning to emerge with evidence of all three canonical resistance-associated SNPs in various hookworm species (27,30,(65)(66)(67)(68).
The mutated models again showed that residue E198 is likely to be a key player in BZ binding in b-tubulins. Similar to previous work on the strongyle parasite H. contortus, in Anc. duodenale, while the F167Y mutation does not directly affect the interaction of BZ, with the bond between BZ and E198 still being seen, it may have an indirect effect (19). We saw a self-binding interaction between the mutated F167Y residue and Q134. The interaction between these two residues has been previously predicted to affect the shape of the entry space of the binding pocket and therefore may affect BZ's ability to enter the binding region, thereby leading to resistance (19). This same self-binding interaction was not seen in the Tr. trichiura 167Y model, nor was it seen in Ascaris b-tubulin models in a previous study (20). This, alongside the evidence of the F167Y SNP in BZ-susceptible populations of Ascaris and the findings of a reduction on the proportion of F167Y SNPs post-treatment in Tr. trichiura, could suggest that the entry blocking interaction between Q134 and F167Y may be a strongyle-specific trait (27).
In both the Anc. duodenale E198A and Tr. trichiura E197A models, we saw a complete lack of binding between the BZ and this key residue. There was also a lack of any bonds of comparable strength in these models, suggesting that the loss of this key residue leads to a less stable binding conformation of BZs. The biggest difference between Anc. duodenale and Tr. trichiura was seen with the F200Y and F199Y SNPs, respectively. In Anc. duodenale, we saw the loss of the key binding interaction with E198, and a self-binding interaction was seen between the mutated F200Y residue and E198. In Tr. trichiura, the binding between E197 and F199Y was also seen, although strong ionic bonds were also formed between BZ and E198. The interaction between the mutated F200Y residue and E198 was similarly seen in other parasites, and it was predicted that it could lead to the destabilization of BZ in the binding pocket and therefore cause resistance (19,20).

CONCLUSION
The results presented here show that BZs act via the same mechanisms in all species tested and once again highlights E198 as a key residue for BZ binding. This analysis includes a wide range of helminths designated as STHs of human importance and closely related parasites of animals. The major difference between the groups of helminths is in the number of btubulin isotypes found within the genome rather than the diversity of these isotypes.
We can postulate from this that the differences seen in the development of resistance are not due to the structural changes in b-tubulin proteins themselves but may be the result of the differences in the b-tubulin isotype abundance and expression. However, it must be noted that the finding of the Q134-F167Y interaction in Anc. duodenale together with evidence of this F167Y SNP in susceptible Ascaris populations and its reduction in frequency after BZ treatment in Tr. trichiura may suggest that mutations in F167Y have more of an impact in strongyle parasite than other STHs.
There is still a need to be vigilant for the development of resistance in ascarid species and other nematodes of veterinary and medical importance. As we are beginning to see reports of treatment failure in ascarids and STHs, genetic investigation will aid in determining the underlying causes in this group of nematodes that have, until recently, been the only group of parasitic nematodes with no significant evidence of BZ resistance.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.