Variation in anthelmintic responses are driven by genetic differences among diverse 1 C. elegans wild strains 2

38 Treatment of parasitic nematode infections in humans and livestock relies on a small arsenal of 39 anthelmintic drugs that have historically reduced parasite burdens. However, anthelmintic 40 resistance (AR) is increasing, and little is known about the molecular and genetic causes of 41 resistance for most drugs. The free-living roundworm Caenorhabditis elegans has proven to be a 42 tractable model to understand AR, where studies have led to the identification of molecular targets 43 of all major anthelmintic drug classes. Here, we used genetically diverse C. elegans strains to 44 perform dose-response analyses across 26 anthelmintic drugs that represent the three major 45 anthelmintic drug classes (benzimidazoles, macrocyclic lactones, and nicotinic acetylcholine 46 receptor agonists) in addition to seven other anthelmintic classes. First, we found that C. elegans 47 strains displayed significant variation in anthelmintic responses across drug classes. Dose-48 response trends within a drug class showed that the C. elegans strains elicited similar responses 49 within the benzimidazoles but variable responses in the macrocyclic lactones and nicotinic 50 acetylcholine receptor agonists. Next, we compared the effective concentration estimates to 51 induce a 10% maximal response (EC 10 ) and slope estimates of each dose-response curve of each 52 strain to the reference strain, N2, which enabled the identification of anthelmintics with population-53 wide differences to understand how genetics contribute to AR. Because genetically diverse strains 54 displayed differential susceptibilities within and across anthelmintics, we show that C. elegans is 55 a useful model for screening potential nematicides. Third, we quantified the heritability of 56 responses to each anthelmintic and observed a significant correlation between exposure closest 57 to the EC 10 and the exposure that exhibited the most heritable responses. Heritable genetic 58 variation can be explained by strain-specific anthelmintic responses within and across drug 59 classes. These results suggest drugs to prioritize in genome-wide association studies, which will 60 enable the identification of AR genes.


ABSTRACT
Treatment of parasitic nematode infections in humans and livestock relies on a small arsenal of anthelmintic drugs that have historically reduced parasite burdens.However, anthelmintic resistance (AR) is increasing, and little is known about the molecular and genetic causes of resistance for most drugs.The free-living roundworm Caenorhabditis elegans has proven to be a tractable model to understand AR, where studies have led to the identification of molecular targets of all major anthelmintic drug classes.Here, we used genetically diverse C. elegans strains to perform dose-response analyses across 26 anthelmintic drugs that represent the three major anthelmintic drug classes (benzimidazoles, macrocyclic lactones, and nicotinic acetylcholine receptor agonists) in addition to seven other anthelmintic classes.First, we found that C. elegans strains displayed significant variation in anthelmintic responses across drug classes.Doseresponse trends within a drug class showed that the C. elegans strains elicited similar responses within the benzimidazoles but variable responses in the macrocyclic lactones and nicotinic acetylcholine receptor agonists.Next, we compared the effective concentration estimates to induce a 10% maximal response (EC10) and slope estimates of each dose-response curve of each strain to the reference strain, N2, which enabled the identification of anthelmintics with populationwide differences to understand how genetics contribute to AR.Because genetically diverse strains displayed differential susceptibilities within and across anthelmintics, we show that C. elegans is a useful model for screening potential nematicides.Third, we quantified the heritability of responses to each anthelmintic and observed a significant correlation between exposure closest to the EC10 and the exposure that exhibited the most heritable responses.Heritable genetic variation can be explained by strain-specific anthelmintic responses within and across drug classes.These results suggest drugs to prioritize in genome-wide association studies, which will enable the identification of AR genes.

AUTHOR SUMMARY
Parasitic nematodes infect most animal species and significantly impact human and animal health.Control of parasitic nematodes in host species relies on a limited collection of anthelmintic drugs.However, anthelmintic resistance is widespread, which threatens our ability to control parasitic nematode populations.Here, we used the non-parasitic roundworm Caenorhabditis elegans as a model to study anthelmintic resistance across 26 anthelmintics that span ten drug classes.We leveraged the genetic diversity of C. elegans to quantify anthelmintic responses across a range of doses, estimate dose-response curves, fit strain-specific model parameters, and calculate the contributions of genetics to these parameters.We found that genetic variation within a species plays a considerable role in anthelmintic responses within and across drug classes.Our results emphasize how the incorporation of genetically diverse C. elegans strains is necessary to understand anthelmintic response variation found in natural populations.These results highlight drugs to prioritize in future mapping studies to identify genes involved in anthelmintic resistance.

INTRODUCTION
Parasitic nematodes are incredibly diverse and infect most animal and plant species [1,2].
Treatments rely on a limited arsenal of anthelmintic drugs with the same drug classes used across most parasite species [3].The three major anthelmintic classes are benzimidazoles (BZs), macrocyclic lactones (MLs), and nicotinic acetylcholine receptor (nAChR) agonists.Over-reliance and inappropriate use of these anthelmintics have placed strong selective pressures on parasites and caused the evolution of anthelmintic resistance (AR) in every drug class [3].Over time, AR causes the drug to become ineffective [4].In many cases, AR is highly heritable, which suggests that natural genetic variants play an important role in the evolution of resistant nematodes [5,6].
With the heavy burden parasitic nematodes place on global health, the limited suite of drugs currently available, and AR on the rise, it is imperative that we identify resistance alleles.With this knowledge, we can responsibly apply drugs and implement treatment strategies to reduce our global infection rate and burden of parasitic nematodes.
Most of what we know about mechanisms of resistance comes from studies of a single strain within a single species, the laboratory-adapted strain of Caenorhabditis elegans called N2.
However, a single genetic background, whether in a free-living or parasitic species, cannot capture the enormous diversity present in the entire species, nor can it predict how natural populations of parasitic nematodes will respond to a drug [7].In aggregate, single strains from many species might capture phylum-level variation, which strengthens the opportunity to identify genes involved in mechanisms of resistance.It is difficult to accurately test AR in genetically diverse parasitic nematodes because of multiple factors: a lack of access to relevant life cycle stages, lack of global sample collections, host-dependent and cost-intensive laboratory life cycles, complex or non-existent in vitro culture systems, and a limited molecular toolkit [7].
With its ease of growth, genetic tractability, and ample molecular toolkit, the roundworm C. elegans is our most useful model to study AR.To date, C. elegans has contributed to the identification and characterization of mechanisms of resistance of all major anthelmintic drug classes [8][9][10][11][12][13][14].Additionally, the natural genetic variation across the C. elegans species is accessible and continuously archived in the C. elegans Natural Diversity Resource (CeNDR), which has facilitated the characterization of natural responses to anthelmintic drugs [11,15,16].
Whole-genome sequence data and identified AR variants are available for CeNDR strains and can be queried to identify orthologous genes between C. elegans and parasites [8].Lastly, because of the tractability of C. elegans and established high-throughput assays (HTA), we can measure C. elegans responses to any soluble compound [17].Thus, the genetic diversity of C. elegans can enable the discovery of the molecular targets of anthelmintics, which are likely to translate across parasitic nematode species.
Here, we performed dose-response analyses that used 26 anthelmintics across six genetically diverse C. elegans strains to identify how growth rate was affected.The anthelmintics used in this study represent the three major anthelmintic drug classes (BZs, MLs, and nAChR agonists) in addition to seven other classes of anthelmintics (nicotinic acetylcholine receptor [nAChR] antagonists, a pore-forming crystal protein, a cyclicoctadepsipeptide, diethylcarbamazine, piperazine, a salicylanilide, and schistosomicides).We measured nematode development after drug exposure for six genetically diverse C. elegans strains using an established high-throughput phenotyping assay [17].We assayed the strains with high levels of replication, collecting a total of 48,343 replicate anthelmintic responses across genetically diverse C. elegans strains, a throughput not possible using parasites.We used phenotypic responses to each anthelmintic to estimate dose-response curves, fit strain-specific model parameters, and calculate the contributions of genetics to these anthelmintic responses.Our results emphasize how the incorporation of natural genetic variation is necessary to quantify drug responses and identify the range of drug susceptibilities in natural populations.Importantly, studies focusing on genetic variation increase the likelihood of identifying orthologous genes between C. elegans and parasites of interest and, in turn, discover mechanisms of resistance shared across species [8].

RESULTS AND DISCUSSION
High-throughput assays across six wild strains facilitated dose-response assessments of

anthelmintic drugs
Dose-response assessments were performed using a microscopy-based high-throughput phenotyping assay for developmental delay in response to 26 anthelmintics (Fig 1).Anthelmintics assayed were five different BZs, seven MLs, three nAChR agonists, three nAChR antagonists, one pore-forming crystal toxin, one cyclooctadepsipeptide, one diethylcarbamazine, one pyrazinoisoquinoline, one salicylanilide, and three schistosomicides (Table 1).Six genetically diverse C. elegans strains were exposed to each anthelmintic in high replication.After measuring nematode responses, phenotypic data were cleaned and processed (see Methods).Next, doseresponse curves were estimated for each anthelmintic to describe how genetic variation contributed to differences in anthelmintic resistance among strains.Differences in responses were measured by the change in developmental rate, as measured by animal length.Nematodes grow longer over time, and anthelmintics have been shown to slow this development [9,10,[18][19][20][21].
Therefore, shorter animals after drug exposure demonstrated that the anthelmintic had a detrimental effect on development.
Figure 1.The high-throughput phenotyping assay allows for rapid dose-response assessments across genetically diverse C. elegans strains.A) strains were passaged for three generations to reduce generational effects of starvation.B) strains were bleach synchronized to collect embryos and then hatched and arrested at the L1 larval stage.C) Animals were fed and exposed to anthelmintics.D) After 48 hours of growth, animals were imaged to collect phenotypic measurements.E) Data were cleaned, and dose-response analysis was performed.Detailed descriptions of all steps can be found in Methods.Created with BioRender.com.Modified from a previous version [22].Table 1.Drug class, subclass, and drugs used in this study.
A four-parameter log-logistic dose-response curve was modeled for each of the 26 anthelmintics, where normalized median animal length was used as the metric for a phenotypic response (see Methods).For each strain-specific dose-response model, slope (b) and effective concentration (e) were estimated with strain as a covariate (S1 Table, S2 Table ).EC10 estimates were found to be more heritable than half maximal effective concentration (EC50) estimates and were therefore used throughout our analyses (S1 Fig).Dose-response relationships described how different strains were affected at varying levels of anthelmintic exposure providing insights into how genetic differences impact anthelmintic susceptibility.
To test for differences in AR among the strains, we looked for differences in the strainspecific dose-response model parameters.We found differences in EC10 values for 22 anthelmintics (S1 Table ).Next, we focused on EC10 comparisons between the reference strain N2 and all other strains (S3 Table ).In total, we observed 44 instances across 22 compounds where at least one strain was significantly more resistant or sensitive than the reference strain N2 using EC10 as a proxy (Student's t-test, Bonferroni correction; padj < 0.05).Because most studies in C. elegans AR have been conducted using the laboratory reference strain, N2, or mutant strains in the N2 genetic background, these results emphasize the importance of using genetically diverse individuals to understand drug responses.Furthermore, the observed frequency of strains with significantly greater anthelmintic sensitivity than the N2 strain was different than what is expected under the null expectation (see Methods; Fisher exact test; p < 0.05), which suggests that diverse C. elegans strains are not equally likely to be susceptible or resistant with respect to the commonly used N2 reference strain.The strain MY16 displayed the most sensitivity and resistance compared to the N2 strain, making up 31% and 63% of the cases, respectively.
Out of the 130 strain-specific slope comparisons with respect to the N2 strain, we observed 92 instances across the 26 compounds where a strain had a significantly different slope than the N2 strain (S4 Table ).However, slope estimate comparisons between the N2 strain and every other strain only describe part of the breadth of C. elegans anthelmintic responses.For this reason, we compared all slopes in a pairwise fashion.Out of the 390 total strain-specific slope comparisons, we observed 275 pairwise instances across the 26 compounds where one strain had a significantly different slope than another strain (S4 Table ).The variation in strain-specific slope comparisons further supports how the incorporation of genetic diversity is necessary to identify anthelmintic responses within a species.Here, we reinforce what previous studies have shown, that C. elegans is a powerful model for assessing the impact of genetic differences on phenotypic variation [23].

Variation in response to BZs is driven by genetic differences among naturally diverse strains
Although BZs are essential in human and veterinary medicine, resistance to this drug class is prominent and common in natural parasite populations [24,25].Historically, the mechanisms of nematode resistance to BZs were thought to have been limited to variants in the drug target betatubulin [26][27][28][29].However, genetic differences in beta-tubulin genes do not explain all intraspecific and interspecific variations observed in BZs efficacy [30] or in responses to different BZs derivatives [31,32].Genome-wide association studies (GWAS) of responses to albendazole, a widely used BZ, found quantitative trait loci (QTL) that do not overlap with beta-tubulin genes, suggesting that additional genes are involved in albendazole resistance [18].Additionally, previous work, which included genetically diverse strains of C. elegans and Caenorhabditis briggsae, a closely related selfing species, found that conserved and drugspecific loci contribute to the effects of BZs (albendazole, fenbendazole, mebendazole, and thiabendazole) [33].Because of evidence that additional genes beyond beta-tubulin genes are involved in BZs resistance, we have yet to fully understand the mechanisms of BZs resistance.
We assessed how natural variation contributes to phenotypic responses across five clinically relevant BZs (albendazole, benomyl, fenbendazole, mebendazole, and thiabendazole) that are widely used in human and veterinary medicine.The panel of six genetically divergent C. elegans wild strains was exposed to increasing concentrations of the five BZs (S5 Table , Fig    2).The strain MY16 displayed resistance in all five BZ dose-response curves, where the EC10 for MY16 was significantly higher than EC10 estimates from all the other strains in every BZ (Fig 2).
The MY16 strain has a non-synonymous variant in the beta-tubulin gene ben-1, causing an amino acid change (A185P) [34] and a presumptive reduction in ben-1 function.The other five strains do not have variants known to reduce ben-1 function.
Benzimidazole strain-specific slope (b) estimates for each dose-response model varied but followed similar trends compared to EC10 estimates (Fig 3A, Fig 3B).These results suggest that the genetic differences among C. elegans strains mediate differential susceptibility across BZs.To quantify the degree of phenotypic variation attributable to segregating genetic differences among strains, we first estimated broad-sense heritability (H 2 ) and narrow-sense heritability (h 2 ) of the phenotypic response for each dose of every BZ (see Methods) (Fig 3C ).For example, we observed that H 2 ranged from 0 in 1 µM albendazole to 0.87 in 51.54 µM albendazole, and h 2 ranged from 0 in 1 µM albendazole to 0.73 in 51.54 µM albendazole.This heritable response indicated that genetic differences among the six strains underlie the variation in albendazole responses.Importantly, all five BZs had highly heritable responses, which indicates that the genetic diversity of C. elegans can be used to identify additional molecular mechanisms of BZs resistance beyond the ben-1 beta-tubulin gene.Because ben-1 is not the only gene involved in BZs responses [18,33,35], we removed the MY16 strain from our analyses to observe how smaller genetic effects play a role in BZs responses in the other five strains (S2 Fig) .After removing MY16, we observed that the strain CX11314 displayed the greatest resistance among the remaining five strains in all BZs except thiabendazole.The strain N2 displayed the greatest resistance in thiabendazole after MY16.Also, the strains CB4856 and JU775, previously described as sensitive to BZs [34,36], displayed sensitivity across BZs and significant variability in thiabendazole, where the JU775 strain was more sensitive than the CB4856 strain (S2 Fig) .Even after removing MY16, we found that responses for thiabendazole were highly heritable, although moderately heritable responses were observed for albendazole.Benomyl, fenbendazole, and mebendazole had reduced heritability.
These results support the previous findings that ben-1 is not the only gene involved in BZs resistance and that diverse C. elegans strains vary across a spectrum of BZ responses [18,33,35].
The strain MY16 is a striking example of how well natural BZ resistance alleles can protect nematodes from BZ treatment.In the context of natural parasitic nematode populations, it is easy to imagine how such beneficial alleles could spread rapidly and further exacerbate parasitic burdens.

Small variations in MLs dosage can significantly alter drug effectiveness among naturally diverse strains
The MLs comprise avermectins and milbemycins and are an essential class of anthelmintics because of our high dependence on them to control nematode parasites in livestock, companion animals, and humans [37].Previous genetic screens performed in the C. elegans laboratory-adapted reference strain, N2, identified three genes that encode glutamategated chloride (GluCl) channel subunits (glc-1, avr-14, and avr-15) that are targeted by MLs [38,39].Studies of abamectin have found additional loci involved in resistance [36].By contrast, ML resistant parasitic nematode isolates do not have mutations in genes that encode GluCl channel subunits, suggesting additional mechanisms of resistance to MLs exist [40,41].
Quantitative genetic mappings in free-living and parasitic nematode species have identified genomic regions that confer drug resistance [18,33,36,[41][42][43][44].Altogether, mutations in GluCl channel genes have modest effects on some ML responses and do not explain the full spectrum of AR in this class.Other genes must be investigated to understand ML mechanisms of resistance.
Here, we assessed how natural variation contributes to phenotypic responses across seven MLs composed of five avermectins (abamectin, doramectin, eprinomectin, ivermectin, and selamectin) and two milbemycins (milbemycin and moxidectin) (S5 Ivermectin did not have significantly different EC10 results among the six strains, suggesting that natural variation in these strains does not affect ivermectin resistance.Moxidectin had undefined EC10 (estimates greater than the maximum exposure) and slope estimates, suggesting higher doses are needed to measure phenotypic responses across strains.Taken together, these results suggest that the genetic differences among C. elegans strains mediate differential susceptibilities across the majority of MLs.
To quantify the degree of phenotypic variation attributable to segregating genetic differences among strains, we estimated the H 2 and h 2 of the phenotypic responses in all MLs (Fig 5C).We observed that H 2 ranged from 0.02 in 0.00533 µM ivermectin to 0.87 in 0.27 µM milbemycin, and h 2 ranged from 0.01 in 0.00105 µM doramectin to 0.73 in 0.27 µM milbemycin.
Heritability for moxidectin could not be calculated because modeling produced undefined EC10 and slope estimates (Fig 4B).In this strain set, we found milbemycin had the highest heritability estimate, whereas ivermectin and selamectin had the lowest heritability estimates of the MLs, indicating genetic variants between the six strains are involved in milbemycin response.
Another important factor in AR is suboptimal dosing of anthelmintics.Misdosing can cause variability in how the drug reaches targeted nematodes and causes an insufficient anthelmintic dose, which allows parasitic nematode populations to develop AR [45].Although error-prone dosing methods can impact AR across all drug classes, it may be particularly important for MLs because small changes can cause vastly different anthelmintic responses.Here, we showed that small changes in MLs doses can significantly vary in effectiveness because of the steep response curves (Fig 4).Additionally, the bioavailability of an anthelmintic and the length of exposure time also play a role in the dosage required to eliminate parasitic nematodes [46].Correct dosing and appropriate bioavailability are critical in all anthelmintic treatments, but this point is even more striking in the MLs where the effective dose range is small.

Natural genetic variation across C. elegans strains explains nAChR agonists responses
The nematode nAChRs in muscle cells are the target of the cholinergic agonists (e.g., levamisole, pyrantel, and morantel [47]).These nAChR agonists cause ligand-gated ion channels to open, producing prolonged muscle contraction and spastic paralysis in nematodes [47].
Here, we assessed how natural variation contributes to phenotypic responses across three nAChR agonists composed of tetrahydropyrimidines (morantel and pyrantel) and an imidazothiazole (levamisole) (S5 estimates for each dose-response model varied but followed near identical trends compared to EC10 estimates, indicating genetic variation is responsible for the observed response variation (Fig 7).The strain DL238 had the highest EC10 in levamisole.We found that the strain CB4856 had the highest EC10 for both tetrahydropyrimidines, whereas the CX11314 and MY16 strains had the lowest EC10 and are thus the most susceptible.Variable patterns across strains within the same drug class suggest nAChR agonists might be acting on different genetic targets.
EC10 and strain-specific slope (b) estimates suggested that the genetic differences among C. elegans strains mediate differential susceptibility across nAChR agonists.To quantify the degree of phenotypic variation attributable to segregating genetic differences among strains, we estimated the H 2 and h 2 of the phenotypic response for each dose of every nAChR agonist (Fig

7C)
. We observed that H 2 ranged from 0.06 in 300 µM pyrantel to 0.52 in 282.5 µM morantel, and h 2 ranged from 0.028 in 300 µM pyrantel to 0.34 in 82.5 µM morantel.We found morantel to elicit the most heritable response, whereas pyrantel elicited the lowest heritable response of the nAChR agonists.Variable heritability in the tetrahydropyrimidines indicates that nAChR agonists may be acting on different genetic targets.Even if drugs have similar trends in EC10 and slope, heritability might help identify drugs where phenotypic variance in response to anthelmintic treatment is attributable to genetic differences.Additionally, although we have several genetic targets identified in C. elegans, it is unclear whether nAChR gene families remain highly conserved among nematode species or whether different species-specific functions can be exploited as potential targets for the control of particular parasites [54].In the past few decades, Cry5B and the nAChR agonist (Levamisole) have been used in combination therapy as strains resistant to nAChR agonists were susceptible to Cry5B [55].Here, we find that the CB4856, DL238, and MY16 strains were sensitive to Cry5B, whereas the We observed different patterns of susceptibility (strain rank order) between levamisole and Cry5B, indicating that this combined therapy could be an effective drug combination.Another promising combination therapy is derquantel and abamectin [56,57].Derquantel and abamectin have been used in combination to treat multi-drug resistant Haemonchus contortus [58,59].However, studies have found monepantel to be more effective than the combined derquantel and abamectin treatment, although monepantel resistance is also prevalent, further exaggerating resistance issues in H. contortus [56,60,61].Here, we observed that the strain MY16 was most sensitive to abamectin, whereas the CB4856 and DL238 strains were most resistant.In derquantel, we found that the strains JU775 and MY16 were the most sensitive.Comparatively, the strain JU775 showed significant sensitivity to the monepantel drugs (monepantel sulfone, monepantel sulfide) Here, patterns of susceptibility and resistance indicated that combination therapy composed of derquantel and abamectin would be a more effective treatment than monepantel alone.Lastly, emodepside has also been commercialized and approved for anthelmintic treatments in companion animals in combination with praziquantel [62].Here, we find that emodepside had heritable responses (S7 Promising combination therapies can be tested using larger strain sets.Altogether, doseresponse assessments in C. elegans provide a useful platform to assess hypothesized effectiveness of drugs that can be used together in combination therapies.Figure 8. Variation in EC10 and dose-response slope estimates can be explained by genetic variation across strains.A) The relative potency of each anthelmintic for each strain compared to the N2 strain is shown.Solid points denote strains with significantly different relative resistance to that anthelmintic compared to the N2 strain (Student's t-test and subsequent Bonferroni correction with a padj < 0.05).Faded points denote strains not significantly different than the N2 strain.Asterisks denote strains with normalized estimates greater than ten compared to the N2 strain.See Supplementary Figure 14 for the relative potency of all strains in each anthelmintic.Anthelmintic drugs with undefined EC10 estimates (estimates greater than the maximum dose to which animals were exposed) are not shown.B) For each anthelmintic, the relative steepness of the dose-response slope inferred for that strain compared to the N2 strain is shown.Solid points denote strains with significantly different dose-response slopes for that anthelmintic compared to the N2 strain (Student's t-test and subsequent Bonferroni correction with a padj < 0.05).Faded points denote strains without significantly different slopes than the N2 strain.Asterisks represent strains with slope estimates greater than 20 compared to the N2 strain.See Supplementary Figure 15 for the slope estimates of all strains in each anthelmintic.Anthelmintic drugs with undefined slope estimates are not shown.The broad class to which each anthelmintic belongs is denoted by the strip label for each facet.

Dose-response assessment across genetically diverse stains identifies five anthelmintics for which C. elegans had little to no phenotypic responses
Because C. elegans is an inexpensive and highly tractable model, we could quickly assess which drugs we should continue to study in this model and drugs that will likely not provide useful results.Dose-response curves for the schistosomicides (niridazole, oxamniquine, and With minimal to no response for oxamniquine, praziquantel, and piperazine, EC10 and slope estimates could not be calculated.The schistosomicides have previously been shown to have no activity against nematodes [63], and thus it was not surprising that we observed little response in C. elegans.Diethylcarbamazine contains a piperazine ring that is essential for the activity of the drug and is the treatment of choice for lymphatic filariasis and loiasis [64].Because piperazine did not cause a response in C. elegans, it is not surprising that diethylcarbamazine did not as well.Although little to no responses (i.e., nematode growth defects) were observed for five drugs using our assays, it is possible assays measuring different fitness traits could elicit anthelmintic effects.Overall, C. elegans is a useful model for screening potential nematicides.

C. elegans is an invaluable model for understanding anthelmintic drug target identification and characterization
Our assays measured C. elegans growth rates in the presence of anthelmintics across multiple concentrations of each drug, and this study yielded several major findings.First, doseresponse experiments captured the variation in growth rate across six diverse strains within the C. elegans species.Growth trends yielded information that can be used to assess how other nematodes might respond to the tested anthelmintic drugs.Second, the large-scale HTA provided quantitative data with the required statistical power and sample sizes needed to effectively measure anthelmintic responses.Third, we were able to identify which drugs had heritable responses and, of those drugs, which doses were most heritable.The most heritable doses of each drug can be used in downstream GWAS to identify genomic regions correlated with AR [15].
By narrowing genomic intervals, subsequent candidate genes can be identified, edited, and validated for anthelmintic responses to ultimately identify the genetic variants involved in resistance and inform downstream parasitic nematode treatments [33,36].
It is well established that parasitic nematodes are more genetically diverse than C. elegans and infect virtually all animal species, so understanding the role genetic diversity plays in anthelmintic resistance is critical [8].The flow of resistance alleles within and among parasite populations has profound implications for the epidemiology of host infection, disease presentation, and the responses of parasite populations to selection pressures, such as anthelmintic treatment [5,65].Additionally, the development of novel anthelmintics is slow, expensive, and complex, making it critical to correctly apply and monitor the usage of our existing drugs.We suggest that genetically diverse C. elegans strains should be deployed to aid highthroughput anthelmintic screening efforts to identify effective anthelmintics and estimated effective concentrations to use when testing in parasitic nematodes.
The presented data focused on the natural genetic variation in C. elegans and will require additional studies to identify genes responsible for the observed anthelmintic responses.This study summarized anthelmintic responses in naturally diverse C. elegans strains and highlighted drugs to focus on in downstream studies.Ultimately, AR gene variants responsible for observed effects need to be identified and validated in C. elegans and subsequently tested in parasitic nematodes.

C. elegans strain selection and maintenance
Six Caenorhabditis elegans strains (CB4856, CX11314, DL238, JU775, MY16, PD1074) from the C. elegans Natural Diversity Resource (CeNDR) were used in this study [16].Isolation details for the six strains are included in CeNDR.These strains were selected from the CeNDR divergent strain set, where the strain PD1074 is referred to by its isotype name, N2.In CeNDR, strains that are >99.97%genetically identical are grouped into isotypes, PD1074 and N2 are nearly genetically identical, therefore, we chose to label PD1074 as N2 to illustrate the response of the canonical laboratory strain N2 [66].Before measuring anthelmintic responses, animals were maintained at 20ºC on 6 cm plates with modified nematode growth medium (NGMA), which contains 1% agar and 0.7% agarose to prevent animals from burrowing.The NGMA plates were seeded with the Escherichia coli strain OP50 as a nematode food source.All strains were grown for three generations without starvation on the NGMA plates before anthelmintic exposure to reduce the transgenerational effects of starvation stress.The specific growth conditions for nematodes used in the high-throughput anthelmintic response assays are described below.

Nematode food preparation
Detailed nematode food preparation steps were followed as previously described [22].
One batch of HB101 E. coli was used as a nematode food source for all assays.Briefly, a frozen stock of HB101 E. coli was used to inoculate and grow a one-liter culture at an OD600 value of 0.001.A total of 14 cultures containing one liter of pre-warmed 1x Horvitz Super Broth (HSB) and an OD600 inoculum grew for 15 hours at 37ºC until cultures were in the late log growth phase.After 15 hours, flasks were removed from the incubator and transferred to 4ºC to arrest growth.Cultures went through three rounds of centrifugation, where the supernatant was removed, and the bacterial cells were pelleted.Bacterial cells were washed and resuspended in K medium.The OD600 value of the bacterial suspension was measured and diluted to a final concentration of OD 600 100 with K medium, aliquoted to 15 ml conicals, and stored at -80ºC for use in the anthelmintic dose-response assays.

Anthelmintic stock preparation
All 26 anthelmintic stock solutions were prepared using either dimethyl sulfoxide (DMSO) or water, depending on the anthelmintic's solubility.Sources, catalog numbers, stock concentrations, and preparation for each anthelmintic are provided (S5 Table ).Anthelmintic stock solutions were prepared, aliquoted, and stored at -20ºC for use in the dose-response assays.

High-throughput anthelmintic dose-response assay
For each assay, populations of each strain were amplified and bleach-synchronized in triplicate.The bleach synchronization was replicated to control for variation in embryo survival and subsequent effects on developmental rates that could be attributed to bleach effects.After bleach synchronization, approximately 30 embryos were dispensed into the wells of a 96-well microplate in 50 µL of K medium.Strains were randomly assigned to columns of the 96-well microplates to vary strain column assignments across the replicate bleaches.Each strain was present in duplicate on each plate.Four replicate 96-well microplates within each of the three bleach replicates for each anthelmintic and control condition tested in the assay were prepared, labeled, and sealed using gas-permeable sealing films (Fisher Cat # 14-222-043).Plates were C. elegans animals from images from the HTA [67].The custom CellProfiler pipeline generates animal measurements by using four worm models: three worm models tailored to capture animals at the L4 larval stage, in the L2 and L3 larval stages, and the L1 larval stage, respectively, as well as a "multi-drug high dose" (MDHD) model, to capture animals with more abnormal body sizes caused by extreme anthelmintic responses.Worm model estimates and custom CellProfiler pipelines were written using the WormToolbox in the GUI-based instance of CellProfiler [23].Next, a custom R package, easyXpress (Version 1.0), was then used to process animal measurements output from CellProfiler [17].These measurements comprised our raw dataset.

Data Cleaning
The presented analysis has been modified from previous work [22].All analyses were performed using the R statistical environment (version 4.2.1)unless stated otherwise.The highthroughout anthelmintic dose-response assay produced thousands of images per experimental block; thus, we implemented a systematic approach to assess the quality of animal measurement data in each well.Several steps were implemented to clean the raw image data using metrics indicative of high-quality animal measurements for downstream analysis.
1) Objects with a Worm_Length > 30 pixels, 100 microns, were removed from the CellProfiler data to (A) retain L1 and MDHD-sized animals and (B) remove unwanted particles [68].Using the Worm_Length > 30 pixels threshold to retain small sensitive animals, more small objects, such as debris, were also retained (see Supplementary Information).
2) R/easyXpress [17] was used to filter measurements from worm objects within individual wells with statistical outliers and to parse measurements from multiple worm models down to single measurements for single animals.
3) The data were visualized by drug, drug concentration, assay, strain, and worm model for two purposes.First, to ensure that each drug, by assay, contained control wells that had a mean_wormlength_um between 600 -800 µm, the size of an L4 animal.If the mean_wormlength_um in the control wells was not between the 600 -800 µm range, then that strain and/or assay were removed for the drug (S17 Fig) .This filter ensured the control wells, DMSO or water, primarily contained L4 animals.Assays and drugs that did not meet the control well mean_wormlength_um criteria and were thus subsequently removed were: abamectin (assay A), derquantel (assay H), niridazole (assay H), eprinomectin (assay I), and piperazine (assay I).Second, we wanted to identify drugs that contained a high abundance of MDHD model objects across all assays and drug concentrations.Drugs with an abundance of objects classified by the MDHD model across assays and concentrations likely contain debris.The MDHD model was removed from the following 13 drugs to limit debris and small objects: albendazole, benomyl, Cry5B, diethylcarbamazine, fenbendazole, mebendazole, morantel, niridazole, oxamniquine, piperazine, praziquantel, pyrantel, and thiabendazole.

4)
We then filtered the data to wells containing between three and forty animals, under the null hypothesis that the number of animals is an approximation of the expected number of embryos originally titered into wells (approximately 30).Given that our analysis relied on well median animal length measurements, we excluded wells with less than three animals to reduce sampling error.

5)
We removed statistical outlier measurements within each concentration for each strain for every anthelmintic drug to reduce the likelihood that statistical outliers influence anthelmintic dose-response curve fits.
6) Next, we removed measurements from all doses of each anthelmintic drug that were no longer represented in at least 80% of the independent assays because of previous data filtering steps.
7) Finally, we normalized the data by (1) regressing variation attributable to assay and technical replicate effects and (2) normalizing these extracted residual values to the average control phenotype.For each anthelmintic drug, we estimated a linear model using the raw phenotype measurement as the response variable, and both assay and technical replicate identity as explanatory variables following the formula median_wormlength_um ~ Metadata_Experiment + bleach using the lm() function in base R. We used the residuals from the linear model to remove the effect of assay and bleach from the raw phenotypes.Next, for each drug, we calculated the mean of residual values in control conditions for each strain in each assay and bleach.Finally, for each drug, strain, assay, and bleach, we subtracted the appropriate mean control values from the model residuals to arrive at our normalized length measurements, which were used in all downstream statistical analyses.These normalized length measurements have the helpful property of being centered on zero in control conditions for each strain and, therefore, control for natural differences in the length of the strains.

Small object removal and data cleaning
In previous analyses, we used Worm_Length > 50 (165 microns) to filter out small objects from data before performing cleaning steps [22].For the anthelmintics, we saw that when applying this filter, high dose concentrations for 12 of the 26 anthelmintics were filtered and removed.
Additionally, the anthelmintic selamectin was entirely removed from the dataset (S18 Fig).
Although a Worm_Length > 50 filtered debris from image data, it also filtered small drug-affected nematodes, which were abundant in this study.To ensure that we captured small drug-affected nematodes across anthelmintics and minimized the amount of retained debris, we altered the animal length threshold to Worm_Length > 30 (100 microns).The threshold Worm_Length > 30 was previously recorded as the smallest animal length of L1 animals after an hour of feeding [68].
To confirm that we were retaining animal objects, we (1) retained the MDHD model for drugs that had small animals present at high dose concentrations (see Methods and Data Cleaning) and ( 2) observed high dose well images to ensure the MDHD model was identifying nematodes.

Figure 2 .
Figure 2. Dose-response curves for benzimidazoles (BZs).Normalized animal lengths (y-axis) are plotted for each strain as a function of the dose of benzimidazole supplied in the highthroughput phenotyping assay (x-axis); Albendazole, Benomyl, Fenbendazole, Mebendazole, and Thiabendazole.Strains are denoted by color.Lines extending vertically from points represent the standard deviation from the mean response.Statistical normalization of animal lengths is described in Methods.

Figure 3 .
Figure 3. Variation in benzimidazole (BZ) EC10 dose-response and slope estimates can be explained by genetic variation across strains.A) Strain-specific EC10 estimates (e) for each benzimidazole are displayed for each strain.Standard errors for each strain-and anthelminticspecific EC10 estimates are shown.B) Strain-specific slope estimates (b) for each benzimidazole are displayed for each strain.Standard errors for each strain-and anthelmintic-specific slope estimate are indicated by the line extending vertically from each point.C) The broad-sense (xaxis) and narrow-sense heritability (y-axis) of normalized animal length measurements were calculated for each concentration of each benzimidazole (Methods; Broad-sense and narrowsense heritability calculations).The color of each cross corresponds to the log-transformed dose for which those calculations were performed.The horizontal line of the cross corresponds to the confidence interval of the broad-sense heritability estimate obtained by bootstrapping, and the vertical line of the cross corresponds to the standard error of the narrow-sense heritability estimate.Because ben-1 is not the only gene involved in BZs responses[18,33,35], we removed

Figure 4 .
Figure 4. Dose-response curves for macrocyclic lactones (MLs).Normalized animal lengths (y-axis) are plotted for each strain as a function of the dose of macrocyclic lactone supplied in the high-throughput phenotyping assay (x-axis).Macrocyclic lactones are organized by A) Avermectins: Abamectin, Doramectin, Eprinomectin, Ivermectin, Selamectin; and B) Milbemycins: Milbemycin and Moxidectin.Strains are denoted by color.Lines extending vertically from points represent the standard deviation from the mean response.Statistical normalization of animal lengths is described in Methods.

Figure 5 .
Figure 5. Variation in macrocyclic lactone (ML) EC10 dose-response and slope estimates can be explained by genetic variation across strains.A) Strain-specific EC10 estimates (e) for each macrocyclic lactone are displayed for each strain.Standard errors for each strain-and anthelmintic-specific EC10 estimates are indicated by the line extending vertically from each point.B) Strain-specific slope estimates (b) for each macrocyclic lactone are displayed for each strain.Standard errors for each strain-and anthelmintic-specific slope estimate are indicated by the line extending vertically from each point.C) The broad-sense (x-axis) and narrow-sense heritability (y-axis) of normalized animal length measurements were calculated for each concentration of each macrocyclic lactone (Methods; Broad-sense and narrow-sense heritability calculations).The color of each cross corresponds to the log-transformed dose for which those calculations were performed.The horizontal line of the cross corresponds to the confidence interval of the broadsense heritability estimate obtained by bootstrapping, and the vertical line of the cross corresponds to the standard error of the narrow-sense heritability estimate.EC10, estimated slope, and heritability could not be calculated for moxidectin and, therefore, not plotted.

Figure 6 .
Figure 6.Dose-response curves for nicotinic acetylcholine receptor (nAChR) agonists.Normalized animal lengths (y-axis) are plotted for each strain as a function of the dose of anthelmintic supplied in the high-throughput phenotyping assay (x-axis); Levamisole, Morantel, and Pyrantel.Strains are denoted by color.Lines extending vertically from points represent the standard deviation from the mean response.Statistical normalization of animal lengths is described in Methods.

Figure 7 .
Figure 7. Variation in nicotinic acetylcholine receptor agonists (nAChR agonists) EC10 dose-response and slope estimates can be explained by genetic variation across strains.A) Strain-specific EC10 estimates (e) for each nicotinic acetylcholine receptor agonist are displayed for each strain.Standard errors for each strain-and anthelmintic-specific EC10 estimates are indicated by the line extending vertically from each point.B) Strain-specific slope estimates (b) for each nicotinic acetylcholine receptor agonist are displayed for each strain.Standard errors for each strain-and anthelmintic-specific slope estimate are indicated by the line extending vertically from each point.C) The broad-sense (x-axis) and narrow-sense heritability (y-axis) of normalized animal length measurements were calculated for each concentration of each nicotinic acetylcholine receptor agonist (Methods; Broad-sense and narrow-sense heritability calculations).The color of each cross corresponds to the log-transformed dose for which those calculations were performed.The horizontal line of the cross corresponds to the confidence interval of the broad-sense heritability estimate obtained by bootstrapping, and the vertical line of the cross corresponds to the standard error of the narrow-sense heritability estimates.
Fig, S16 Fig) among genetically diverse C. elegans, but praziquantel had no heritable responses (S13 Fig).Although we describe anticipated anthelmintic effectiveness in combination therapies, we acknowledge that a limitation of this study is the small number of strains used.Even though we used a genetically divergent strain set, we have captured a fraction of the genetic variation across the C. elegans species.
Fig).Minimal responses were also observed for diethylcarbamazine and piperazine (S6 Fig,S12 placed in humidity chambers to incubate overnight at 20ºC while shaking at 170 rpm (INFORS HT Multitron shaker).The following morning, food was prepared to feed the developmentally arrested first larval stage animals (L1s) using the required number of OD600100 HB101 aliquots (see Nematode food preparation).The aliquots were thawed at room temperature, combined into a single conical tube, and diluted to an OD60030 with K medium.To inhibit further bacterial growth and prevent contamination, 150 µM of Kanamycin was added to the HB101.Working with a single anthelmintic at a time, an aliquot of anthelmintic stock solution thawed at room temperature (see Anthelmintic stock preparation) and was diluted to a working concentration.The anthelmintic working concentration was set to the concentration that would give the highest desired dose when added to the 96-well microplates at 1% of the total well volume.The serial dilution of the anthelmintic working solution was prepared using the same diluent, DMSO or water, used to make the stock solution.The dilution factors ranged from 1.2 to 2.5 depending on the anthelmintic used, but all serial dilutions had eight concentrations, including a 0 µM control.The serial dilution was then added to an aliquot of the OD60030 K medium at a 3% volume/volume ratio.Next, 25 µl of the food and anthelmintic mixture was transferred into the appropriate wells of the 96-well microplates to feed the arrested L1s at a final HB101 concentration of OD60010 and expose L1 larvae to an anthelmintic at one of eight levels of the dilution series.Immediately afterward, the 96-well microplates were sealed using a new gas permeable sealing film, returned to the humidity chambers, and incubated for 48 hours at 20ºC shaking at 170 rpm.The remaining 96-well microplates were fed and exposed to anthelmintics in the same manner.After 48 hours of incubation in the presence of food and anthelmintic, the 96-well microplates were removed from the incubator and treated with 50 mM sodium azide in M9 for 10 minutes to paralyze and straighten nematodes.Images of nematodes in the microplates were immediately captured using a Molecular Devices ImageXpress Nano microscope (Molecular Devices, San Jose, CA) using a 2X objective.The ImageXpress Nano microscope acquires brightfield images using a 4.7 megapixel CMOS camera and stores images in a 16-bit TIFF format.The images were used to quantify the development of nematodes in the presence of anthelmintics as described below (seeData collection and Data cleaning).Data processingCellProfiler software (Version 4.0.3) was used to quantify animal lengths from images collected on the Molecular Devices ImageXpress Nano microscope (Carpenter et al. 2006).A Nextflow pipeline (Version 20.01.0) was written to run command-line instances of CellProfiler in parallel on the Quest High-Performance Computing Cluster (Northwestern University).The CellProfiler workflow can be found at (https://github.com/AndersenLab/cellprofiler-nf).CellProfiler modules and Worm Toolbox were developed to extract morphological features of individual