Molecular evidence of widespread benzimidazole drug resistance in Ancylostoma caninum from domestic dogs throughout the USA and discovery of a novel isotype-1 β-tubulin benzimidazole resistance mutation

Ancylostoma caninum is an important zoonotic gastrointestinal nematode of dogs worldwide and a close relative of human hookworms. We recently reported that racing greyhounds in the USA are infected with A. caninum that are commonly resistant to multiple anthelmintics. Benzimidazole resistance in A. caninum in greyhounds was associated with a high frequency of the canonical F167Y(TTC>TAC) isotype-1 β-tubulin mutation. In this work, we show that benzimidazole resistance is remarkably widespread in A. caninum from domestic dogs across the USA. First, we identified and showed the functional significance of a novel benzimidazole isotype-1 β-tubulin resistance mutation, Q134H(CAA>CAT). Several benzimidazole resistant A. caninum isolates from greyhounds with a low frequency of the F167Y(TTC>TAC) mutation had a high frequency of a Q134H(CAA>CAT) mutation not previously reported from any eukaryotic pathogen in the field. Structural modeling predicted that the Q134 residue is directly involved in benzimidazole drug binding and that the 134H substitution would significantly reduce binding affinity. Introduction of the Q134H substitution into the C. elegans β-tubulin gene ben-1, by CRISPR-Cas9 editing, conferred similar levels of resistance as a ben-1 null allele. Deep amplicon sequencing on A. caninum eggs from 685 hookworm positive pet dog fecal samples revealed that both mutations were widespread across the USA, often at high frequency, with prevalences of 49.7% (overall frequency 54.0%) and 31.1% (overall frequency 16.4%) for F167Y(TTC>TAC) and Q134H(CAA>CAT), respectively. Canonical codon 198 and 200 benzimidazole resistance mutations were absent. The F167Y(TTC>TAC) mutation had a significantly higher prevalence and frequency in Western USA than in other regions, which we hypothesize is due to differences in refugia. This work has important implications for companion animal parasite control and the potential emergence of drug resistance in human hookworms. Author Summary Although increasingly common in livestock, no reports of widespread anthelmintic resistance are confirmed in any companion animal or human gastrointestinal nematode parasite to date. The canine hookworm is a common intestinal zoonotic parasite of dogs with severe clinical impacts in young dogs, and for which control is dependent on regular anthelmintic use. We recently reported multiple anthelmintic drug resistance in A. caninum isolates from greyhounds derived from multiple locations in the USA likely caused by long standing intensive treatment regimens in kennels. In this study, we investigated benzimidazole resistance in A. caninum in pet dogs across the USA. We also identified and showed the functional significance of a novel benzimidazole isotype-1 β-tubulin resistance mutation in A. caninum from greyhounds that has not been previously reported in the field for any organism. We then determined that this novel mutation, as well as a previously characterized resistance mutation, were present, often at high frequency, in many A. caninum populations across the USA. This study reports the first evidence of widespread drug resistance for any parasitic nematode of companion animals and illustrates the power of molecular approaches to rapidly assess anthelmintic resistance in a region.


Abstract
Ancylostoma caninum is an important zoonotic gastrointestinal nematode of dogs worldwide and a close relative of human hookworms. We recently reported that racing greyhounds in the USA are infected with A. caninum that are commonly resistant to multiple anthelmintics. Benzimidazole resistance in A. caninum in greyhounds was associated with a high frequency of the canonical F167Y(TTC>TAC) isotype-1 β-tubulin mutation. In this work, we show that benzimidazole resistance is remarkably widespread in A. caninum from domestic dogs across the USA. First, we identified and showed the functional significance of a novel benzimidazole isotype-1 β-tubulin resistance mutation, Q134H(CAA>CAT). Several benzimidazole resistant A. caninum isolates from greyhounds with a low frequency of the F167Y(TTC>TAC) mutation had a high frequency of a Q134H(CAA>CAT) mutation not previously reported from any eukaryotic pathogen in the field. Structural modeling predicted that the Q134 residue is directly involved in benzimidazole drug binding and that the 134H substitution would significantly reduce binding affinity. Introduction of the Q134H substitution into the C. elegans β-tubulin gene ben-1, by CRISPR-Cas9 editing, conferred similar levels of resistance as a ben-1 null allele. Deep amplicon sequencing on A. caninum eggs from 685 hookworm positive pet dog fecal samples revealed that both mutations were widespread across the USA, often at high frequency, with prevalences of 49.7% (overall frequency 54.0%) and 31.1% (overall frequency 16.4%) for F167Y(TTC>TAC) and Q134H(CAA>CAT), respectively. Canonical codon 198 and 200 benzimidazole resistance mutations were absent.
The F167Y(TTC>TAC) mutation had a significantly higher prevalence and frequency in Western USA than in other regions, which we hypothesize is due to differences in refugia. This work has important implications for companion animal parasite control and the potential emergence of drug resistance in human hookworms.

Introduction
Ancylostoma caninum is one of the most important parasitic gastrointestinal nematodes of dogs worldwide, with severe infections causing anemia, and even death in younger puppies and is also of zoonotic importance [1][2][3][4]. It is also an important model for human hookworm infection [5,6]. Treatment and control are dependent on the routine use of broad-spectrum anthelmintic drugs. Although anthelmintic resistance has been widespread in gastrointestinal nematodes of grazing livestock for many years, it has generally not been considered a problem for gastrointestinal nematodes of domestic dogs [4]. Until recently, the only published reports of anthelmintic resistance in A. caninum were several cases of pyrantel resistance in dogs in Australia [7][8][9]. Historically, there has been little concern about the risk of resistance emerging in canine gastrointestinal parasites or the need for drug stewardship with respect to resistance mitigation.
In spite of the historical lack of concern regarding drug resistance, the management of A. caninum has been problematic in greyhounds in the USA for some time, with high infection intensities and frequent reports of cases refractory to anthelmintic treatments. The high hookworm infection intensities in greyhounds are due, at least in part, to the housing conditions in kennels with large numbers of dogs and outdoor runs that favor hookworm transmission. This situation has led to an overreliance on frequent anthelmintic drug treatments to achieve hookworm control, leading to high drug selection pressure being applied to A. caninum, especially in greyhounds in racing and adoption kennels [10]. A recent study confirmed that resistance to multiple anthelmintic drug classes is now common in hookworms in greyhounds originating from kennels of eight different states of the USA; in vitro assays confirmed that all A. caninum isolates that were exposed to benzimidazoles and macrocyclic lactones were resistant [11].
Until very recently, there has been less concern about anthelmintic resistance in A.
caninum from pet dogs as they are subject to much less intense anthelmintic drug selection pressure than kenneled greyhounds. However, there is strong evidence of an increase in the prevalence of A. caninum in pet dogs: An analysis of over 39 million diagnostic reports in pet dogs reported an overall increase of 47% from 2012 to 2018 in hookworm prevalence from 2.02% in 2012 to 2.96% in 2018 [12]. Additionally, more than 3000 fecal samples collected from 288 off-leash dog parks across the USA in 2019, showed a prevalence of 7.1% for A.
caninum [13]. Concurrent with this increase in prevalence, there appears to be a rise in cases reported by veterinarians of persistent A. caninum infections despite multiple rounds of treatment with different drugs. Further, a recent investigation of three canine hookworm cases in pet dogs that were refractory to treatment, showed compelling evidence of multiple anthelmintic drug resistance to benzimidazoles, macrocyclic lactones, and pyrantel using both in vitro and in vivo phenotypic assays [11,14,15], as well as in vivo efficacy based on a controlled clinical trial [16]. However, in contrast to greyhounds, no work has been done to investigate the prevalence of drug resistance in A. caninum from pet dogs in any country .
Hence, we undertook this work to investigate the prevalence and distribution of anthelmintic resistance in hookworms infecting pet dog populations in the USA.
Our current understanding of the molecular genetics of benzimidazole resistance in the strongylid nematode group enables a molecular approach to be used to investigate the prevalence and distribution of resistance for this particular drug class. Molecular screening for drug resistant parasites offers a number of advantages over both in vitro and in vivo phenotypic assays. This allows for larger-scale studies, is more amenable to standardization, and provides information about the molecular epidemiology of resistance. Recent genome-wide studies have confirmed that the single most important benzimidazole resistance locus in the Haemonchus contortus genome is the isotype-1 β-tubulin gene [17,18] [19][20][21][22][23][24][25]. Further, these mutations have all been shown to be sufficient to confer similar levels of benzimidazole resistance in C. elegans by CRISPR-Cas9 editing [26][27][28]. The isotype-1 β-tubulin F167Y(TTC>TAC) mutation, but not mutations at codons 198 and 200, was found to be present at high frequency in all three confirmed multi-resistant A. caninum isolates reported from pet dogs but not in susceptible isolates [14]. The importance of this mutation was further confirmed by a significant correlation between the egg hatch assay IC95 and the F167Y(TTC>TAC) mutation in the greyhounds [11,14].
Here, we use a molecular genetic approach to investigate benzimidazole resistance in A. caninum isolated from pet dogs across the USA. We first report the discovery of a novel benzimidazole resistance mutation at codon 134 in the isotype-1 β-tubulin gene in A. caninum that has not been previously reported in any parasite or fungal species in the field. Genotype-phenotype relationships from greyhounds suggest that benzimidazole resistance can be accounted for by the presence of either the F167Y(TTC>TAC) or the Q134H(CAA>CAT) mutations. Predictive structural modeling shows that the codon 134Q(CAA) residue is important for drug binding and is disrupted by a Histidine(H) substitution at this position. We further show, by CRISPR-Cas9 editing, that the Q134H(CAA>CAT) mutation is sufficient to confer similar levels of benzimidazole resistance in C. elegans as a null ben-1 deletion mutation as with previously reported canonical mutations at codons 167, 198, and 200. Deep amplicon sequencing showed that both the F167Y(TTC>TAC) and Q134H(CAA>CAT) mutations are widespread in A. caninum from pet dogs across the USA, with the former being present at higher frequencies overall. These results demonstrate benzimidazole resistance is well established in A. caninum from pet dogs across the whole of the USA which has important implications for the management of hookworms in dogs and drug stewardship. In addition, the results provide important information for the development and application of molecular diagnostics for benzimidazole resistance and also provide an important precedent for the emergence of resistance in other hookworm species including those of humans. Finally, this work illustrates the value of combining work in parasites with C. elegans for the discovery and confirmation of novel anthelmintic resistance mutations [29].

Discovery of a novel benzimidazole resistance mutation in nematodes: isotype-1 βtubulin Q134H(CAA>CAT)
An Egg Hatch Assay (EHA) was previously conducted on a set of A. caninum isolates from USA greyhounds, and all showed the presence of benzimidazole resistance [11] ( Fig   1A). The canonical isotype-1 β-tubulin F167Y(TTC>TAC) benzimidazole resistance mutation was present at a high frequency in 13 out of the 16 isolates (82.6%, range: 40.2% -98.1%) but was absent in one isolate and was present at low frequency in the remaining two isolates (5.7%, 5.3%) (Fig 1B). The other canonical isotype-1 β-tubulin benzimidazole resistance mutations at codons 198 and 200 were absent from all 16 A. caninum isolates from greyhounds suggesting the potential role of a novel mutation in conferring benzimidazole resistance phenotypes in these three isolates. Examination of the Amplicon Sequence Variants (ASVs) of the 293 bp isotype-1 β-tubulin fragment previously generated by deep amplicon sequencing revealed a single non-synonymous mutation Q134H(CAA>CAT) in many of the greyhound A. caninum isolates and at a high frequency (90.0%, 94.4%, and 100%) in the three isolates in which the F167Y(TTC>TAC) was absent or at low frequency ( Fig 1C). In order to determine if any additional non-synonymous mutations occurred elsewhere in the isotype-1 β-tubulin gene for the haplotype carrying the Q134H(CAA>CAT) mutation, we PCR-amplified and cloned a 3,467 bp fragment encompassing the near fulllength isotype-1 β-tubulin gene from two of the greyhound isolates with a high frequency of this mutation as well as from the benzimidazole-susceptible isolate (Barrow) as a control.
Sequencing of six clones from the three isolates revealed no additional non-synonymous mutations in the isotype-1 β-tubulin (S1 Fig). Consequently, we hypothesized the Q134H(CAA>CAT) to be a novel mutation responsible for the benzimidazole resistance phenotype in the three isolates with a low frequency of the F167Y(TTC>TAC) mutation.
We undertook predictive structural modeling of the A. caninum isotype-1 β-tubulin polypeptide to assess the potential impact of the Q134H(CAA>CAT) mutation on benzimidazole drug binding (Fig 2). The amino acid positions associated with benzimidazole resistance -134, 167, 198, and 200 -were mapped to the nocodazole binding pocket of βtubulin and all of these sites, including the newly discovered Q134(CAA), are in direct Van der Waals contact with the methyl ester terminus of nocodazole. Clinically relevant molecules from the same drug class, such as albendazole, mebendazole, and fenbendazole, vary only at the opposite end of the molecule. Disruption of these contacts by mutation is predicted to introduce steric clashes (Q134H, F167Y, F200Y) or lead to the loss of favorable interactions with the azole ring (E198L/A/V) (Fig 2).
In order to investigate whether the Q134H(CAA>CAT) mutation was sufficient to confer phenotypic resistance in vivo, we performed CRISPR-Cas9 genome editing to introduce a 134H(CAC) allele into the C. elegans ben-1 gene in the N2 genetic background. Two independently edited strains, ean243 and ean244, were created to control for any off-target effects from genome editing. We measured the responses to albendazole (as well as DMSO as vehicle control) exposure of the two ben-1 134H(CAC) genome-edited strains, as well as the N2 (susceptible control) and ben-1 deletion (resistant control) strains . The 134H(CAC) ben-1 allele showed significantly increased levels of resistance to albendazole relative to the susceptible N2 strain and similar levels of resistance as the ben-1 deletion strain (ean64) (Fig   3). These results show that the Q134H(CAC) allele confers similar levels of resistance in C.
elegans as a ben-1 null allele, the same level as has been previously shown for the other ben-1 resistance mutations, including the F167Y(TTC>TAC) mutation [27,28].

USA.
The F167Y(TTC>TAC) isotype-1 β-tubulin benzimidazole resistance mutation has previously been shown to be common in A. caninum infecting greyhounds from racing and adoption kennels in the USA [11]. In order to investigate the situation in pet dogs, hookworm eggs were collected from individual pet dog fecal samples across the USA previously identified as hookworm positive by microscopy. Genomic DNA preparations were made from the hookworm eggs for 328 individual dogs for which a minimum of 300 hookworm eggs were recovered (range 300 -14,000). In order to also examine dogs with low infection intensities (fewer than 300 eggs recovered), DNA lysates were prepared from 65 hookworm egg pools comprising 357 dogs, containing a minimum of 200 eggs per pool (range 200 -10,000) (S1 Appendix, S1 Table). Deep sequencing was performed on two PCR amplicons generated from these DNA preparations -a 293 bp fragment encompassing codons 134 and 167, and a 340 bp fragment encompassing codons 198 and 200 -and variant calling was used to determine the frequency of resistance mutations at each of these codons (Fig 4).  48.7% -59.4%) and a frequency of >50% was present in 87 of these samples (Fig 4A,B).
There was a significant difference in the overall frequency of the 167Y(TAC) resistance mutation across the regions (Kruskal-Wallis, p = 1e-04). The prevalence and overall frequency of the F167Y(TTC>TAC) mutation were significantly higher in the west than in the other three regions of the USA (75% prevalence and 72.8% overall frequency in positive samples from the west) ( Fig 4D, Table 1). The frequency of this mutation was also significantly higher in large breeds of dogs compared to small and medium-sized breeds based on the data from 280 dogs for which breed information was available (Kruskal-Wallis, p = 0.002) (S2 FigA). Out of the 206 dogs with age-related information, the prevalence of the 167Y(TAC) mutation was significantly higher in puppies when compared to seniors (Fisher's test, p = 0.02), but there were no other differences in the frequency of the mutation among the different age groups (S2 FigB).  (Fig 5A,B).
The prevalence and overall frequency of the Q134H(CAA>CAT) resistance mutation were significantly higher in the west compared to the midwest but not compared to the other two regions (overall comparison, Kruskal-Wallis, p = 0.031) ( Fig 5C,   We also examined previously generated amplicon sequencing data from A. caninum from greyhounds sourced from southern USA adoption and racing kennels for the presence of the Q134H(CAA>CAT) resistance mutation [11]. Variant calling at codon 134 on these data revealed that 88.6% (62/ 70) of the A. caninum samples from the greyhounds carried the novel Q134H(CAA>CAT) resistance mutation. Its overall frequency in these positive samples was 25.8% (95% C.I. 20.1% -31.5%) and it was at a frequency of >50% in eight dogs (Fig 6). The prevalence and infection intensities of A. caninum in pet dogs are highest in the southern USA, and generally lower in the west and the northeast regions [30,31]. This difference is likely a consequence of regional differences in ambient temperatures and humidity that favor the environmental transmission of the parasite [32,33]. An interesting result from our study is that the prevalence and frequency of both the F167Y(TTC>TAC) and Q134H(CAA>CAT) mutations are higher in A. caninum isolates from the western USA than in the other regions. This is particularly interesting since the western region of the USA is known to have the lowest prevalence of A. caninum [13,30]. One possible explanation for the higher resistance prevalence and allele frequency is that there is a lower level of environmental refugia in the west due to the drier climate. Low levels of refugia lead to higher drug selection pressure and is thus recognized as an extremely important factor for the development of anthelmintic resistance in nematodes of livestock. This circumstance would therefore parallel the situation with H. contortus in sheep in drier and/or colder climates [34][35][36].

Implications of widespread benzimidazole resistance in
Whether Mass Drug Administration (MDA) programs will lead to selection for benzimidazole resistance in human hookworm species has been a question for many years.
[ [37][38][39]. Our results in A. caninum have some important implications for this issue. First, the study demonstrates that isotype-1 β-tubulin mutations underlie benzimidazole resistance in a parasite that is very closely related to human hookworm species. Second, the work demonstrates the insidious nature of resistance emergence and spread: Our data conclusively show that benzimidazole resistance has become extremely widespread in canine hookworms in the USA, and this occurred without a single previous report of resistance. Given the widespread geographical distribution and high prevalence of resistance we found, the process of spread must have been occurring for many years. The veterinary and scientific communities only recently began debating the cause of apparent hookworm treatment failures, and this only occurred after resistance to all classes of anthelmintics used in dogs had evolved [14].
This emphasizes the importance of pre-emptive surveillance if resistance is to be identified at an early enough stage for mitigation strategies to be of benefit. Thirdly, the model of intense drug treatments in kenneled greyhounds, providing the initial selection pressure for resistance emergence followed by subsequent migration and spread through the pet dog population, has some potentially important lessons for human hookworm control. Once resistance emerges regionally as a result of local drug selection, subsequent spread can occur undetected until at an advanced stage if active surveillance is not undertaken.
Finally, the results presented here provide a basis for the development and application of molecular diagnostic tests for benzimidazole resistance in canine hookworms. The data suggest that the determination of the frequency of both the F167Y(TTC>TAC) and Q134H(CAA>CAT) mutations is necessary for molecular diagnosis and surveillance of benzimidazole resistance in this parasite (Fig 1). Given the wide range in the frequency of these mutations present in different isolates, diagnostic tests will need to be at least semiquantitative to have predictive value for likely treatment success.

The isotype-1 β-tubulin Q134H(CAA>CAT) is a novel benzimidazole resistance mutation
We previously reported that the canonical F167Y(TTC>TAC) isotype-1 β-tubulin Examination of the ASVs from deep sequencing identified the presence of an additional nonsynonymous mutation, Q134H(CAA>CAT), at varying frequencies in 15 out of the 16 isolates ( Fig 1C). Notably, this novel Q134H(CAA>CAT) mutation was present at a high frequency (≥90%) in all three A. caninum isolates (BH25, DL13, DL14) that were benzimidazole resistant but had a low frequency of the F167Y(TTC>TAC) isotype-1 β-tubulin mutation (<6%) (Fig   1B,C). Sequencing of the near full-length isotype-1 β-tubulin gene revealed the Q134H(CAA>CAT) mutation to be the only non-synonymous mutation on the same haplotype.
Protein-drug interaction modeling shows that the residue at codon 134 is directly in contact with the constant portion of the benzimidazole scaffold. By disrupting the contact with the wildtype glutamine, the Q134H(CAA>CAT) mutation would likely significantly reduce benzimidazole drug binding affinity (Fig 2). CRISPR-Cas9 editing in C. elegans showed that the Q134H(CAA>CAC) mutation is sufficient to confer a similar level of benzimidazole resistance in vivo as the previously described canonical resistance mutations at codons 167, 198, and 200 (Fig 3) [27]. These data indicate that Q134H(CAA>CAT) is a resistance mutation that contributes, along with the F167Y(TTC>TAC) mutation, to benzimidazole resistance in A.
caninum in the USA greyhounds. The only previous evidence to suggest the potential of substitutions at codon 134 of a β-tubulin gene to confer benzimidazole resistance is from laboratory mutagenesis experiments on Aspergillus nidulans in which the Q134L mutant strain was benomyl resistant [40]. However, to our knowledge, no mutations at codon 134 have been previously identified or associated with benzimidazole resistance in natural field populations of any organism. From the extensive work that has been conducted on strongylid parasitic nematode species to date, resistance mutations have only been reported at codons 167, 198, and 200 of the β-tubulin gene, codons clustered around the drug binding site [19][20][21][22]25,41]. It is notable that a substitution at codon 134 is now added to this list of benzimidazole resistance mutations found in parasitic nematode field populations in spite of not having been identified previously in extensive studies in C. elegans or fungi. Work from C. elegans also suggests that additional resistance mutations in the isotype-1 β-tubulin gene and perhaps other loci could also be involved in different parasite species and isolates [26,42,43]. These observations not only highlight the importance of directly studying target organisms of interest in the field but also the value of integrating field research on parasitic nematodes with functional analysis in a model organism, in this case, C. elegans. The dogs from which the samples were collected were categorized into small, medium, and large breeds according to the American Kennel Club (AKC, 2021) and into age groups as per recent American Animal Hospital Association guidelines: puppy: less than 1 year of age, young adult: 1 to 3 years, mature adult: 4 to 6 years, and senior: greater than 6 years old [45].

Fecal floatation and DNA extraction
Eggs from hookworm positive fecal samples classified as "many", "moderate", and "few" were isolated individually. Samples classified as "rare", as well as those classified as "few" but with <1 g of feces available, were combined into pools before isolating the eggs.
Wooden tongue depressors were used to mix the fecal samples with water to produce a fecal slurry, which was then filtered through a bleached grade 40 cheesecloth (GRAINGER, Lake University of Calgary). Three, two, and one clones were sequenced from BH45, BH31, and Barrow isolate, respectively, using primers for Sanger sequencing (S2 Table).

CRISPR-Cas9 editing of the Q134H(CAA>CAT) allele into the C. elegans ben-1 gene
C. elegans were cultured on plates of modified nematode growth media (NGMA), containing 1% agar and 0.7% agarose, seeded with OP50 bacteria [46]. Plates were maintained at 20°C for the duration of this experiment. Before the assay, animals were grown for three generations to reduce the multigenerational effects of starvation [47]. The wild-type laboratory reference strain N2 and the previously created ben-1(ean64) deletion strain (ECA882) [26] were used to verify quantitative sensitivity and resistance, respectively. Two strains containing the Q134H(CAA>CAC) allele were generated in the N2 background using CRISPR-Cas9 genome editing as previously described [26,27]. Briefly, after injection, possibly edited strains underwent multiple generations of confirmation, ensuring that strains were homozygous for the Q134H(CAA>CAC) allele. Two independently edited strains, ben-1(ean243) and ben-1(ean244), were created to control for any potential off-target effects.
A previously described high-throughput fitness assay was used for all measurements of benzimidazole responses [26][27][28]. After growing, animals were treated with 50 mM sodium azide in M9 buffer to straighten the animals for scoring [50]. Each individual was measured using the COPAS BIOSORT large particle flow cytometer (Union Biometrica, Holliston MA). Animal optical density normalized by animal length (median.EXT) was calculated for each individual. Optical density is reduced in susceptible animals after exposure to benzimidazoles and is able to differentiate resistant and susceptible animals [27,28]. Median optical density of populations in each well was used for quantification of albendazole response. Higher regressed optical density values indicate reduced anthelmintic response, and lower values indicate a more severe response to anthelmintic exposure.
For the analysis of optical density, all data processing was performed using the R package easysorter [51] as described previously [26,27]. Contaminated wells and outliers were removed, and the data were regressed using the regress() function in the easysorter package to remove control phenotypes from the experimental phenotypes, leaving residual values for analysis. Statistical comparisons were performed in R (Version 4.1.2). Tukey's Honest Significant Difference (HSD) test was performed using the rstatix package [52] on an ANOVA model with the formula "phenotype~strain" to calculate differences among population phenotypes.
The structural model for the A. caninum isotype-1 β-tubulin was made using AlphaFold2 [53] and was then aligned to Porcine β-tubulin bound to nocodazole (PDB 5CA1) using Pymol Version 2.4 [54]. The canonical as well as the novel isotype-1 β-tubulin codons associated with benzimidazole resistance (134, 167, 198, and 200) were inspected and are in direct van der Waals contact with the constant portion of benzimidazole drugs.

Deep amplicon sequencing and analysis
Deep amplicon sequencing was used to identify and determine the frequency of sequence polymorphisms at codons 134, 167, 198, and 200 of the A. caninum isotype-1 βtubulin as previously described [11,14,25]. Two separate fragments, a 293 bp fragment encompassing codons 134 and 167, and a 340 bp fragment, encompassing codons 198 and 200, were PCR-amplified as previously described [11]. Sample purification, indexing, and library preparation for Illumina deep sequencing were performed as previously described [25].
Libraries were sequenced on the Illumina MiSeq platform using the 2 x 250 bp V2 sequencing chemistry (Illumina Inc., USA). One positive control and eight negative controls comprising a benzimidazole-susceptible A. caninum isolate and molecular-grade water, respectively, were included on each 96-well plate.
Following sequencing, the paired-end reads were analyzed using the previously described bioinformatics pipeline with slight modifications [11,55]. As per the developer's recommendation, the reads from each run were processed separately using the DADA2 pipeline, to reduce the error rates. The following parameters were used in addition to the default settings during the quality filtering step in the DADA2 pipeline: i) based on read quality, the forward and reverse reads were trimmed to a length of 220 bp and 150 bp, respectively, for the fragment containing codons 134 and 167; whereas the forward and reverse reads for the fragment containing codons 198 and 200 were trimmed to 220 bp and 200 bp, respectively; ii) reads shorter than 50 bp or with an expected error of >1 and >2 in the forward and reverse reads, respectively, were removed. The filtered and merged reads from each run were then combined, the chimeric sequences were removed using DADA2 and variant analysis was performed as previously described [11]. Additionally, we excluded all samples with a read depth lower than 1000 reads from further analysis. Results of the variant analysis were visualized in R (Version 4.1.0) using the packages ggplot2 and leaflet [56,57].

Assessment of the repeatability of SNP frequencies determined by amplicon sequencing
In order to test the repeatability of the benzimidazole resistance SNP frequencies determined by amplicon sequencing, two independent PCRs were performed on 50 randomly    -1(ean64) is a previously created deletion strain. ben-1(ean243) and ben-1(ean244) are two independent CRISPR-Cas9 edited Q134H(CAA>CAC) alleles. Each point represents the regressed phenotype of a well containing approximately 50 animals. Data are shown as Tukey box plots with the median as a solid horizontal line, and the top and bottom of the box representing the 75th and 25th quartiles, respectively. The top whisker is extended to the maximum point within the 1.5 interquartile range from the 75th quartile. The bottom whisker is extended to the minimum point within the 1.5 interquartile range from the 25th quartile. Statistical significance is shown as asterisk above each strain (p < 0.001 = ***, p < 0.0001 = ****, Tukey HSD).