Transcriptomic analyses implicate neuronal plasticity and chloride homeostasis in ivermectin resistance and recovery in a parasitic nematode

The antiparasitic drug ivermectin plays an essential role in human and animal health globally. However, ivermectin resistance is widespread in veterinary helminths and there are growing concerns of sub-optimal responses to treatment in related helminths of humans. Despite decades of research, the genetic mechanisms underlying ivermectin resistance are poorly understood in parasitic helminths. This reflects significant uncertainty regarding the mode of action of ivermectin in parasitic helminths, and the genetic complexity of these organisms; parasitic helminths have large, rapidly evolving genomes and differences in evolutionary history and genetic background can confound comparisons between resistant and susceptible populations. We undertook a controlled genetic cross of a multi-drug resistant and a susceptible reference isolate of Haemonchus contortus, an economically important gastrointestinal nematode of sheep, and ivermectin-selected the F2 population for comparison with an untreated F2 control. RNA-seq analyses of male and female adults of all populations identified high transcriptomic differentiation between parental isolates, which was significantly reduced in the F2, allowing differences associated specifically with ivermectin resistance to be identified. In all resistant populations, there was constitutive upregulation of a single gene, HCON_00155390:cky-1, a putative pharyngeal-expressed transcription factor, in a narrow locus on chromosome V previously shown to be under ivermectin selection. In addition, we detected sex-specific differences in gene expression between resistant and susceptible populations, including constitutive upregulation of a P-glycoprotein, HCON_00162780:pgp-11, in resistant males only. After ivermectin selection, we identified differential expression of genes with roles in neuronal function and chloride homeostasis, which is consistent with an adaptive response to ivermectin-induced hyperpolarisation of neuromuscular cells. Overall, we show the utility of a genetic cross to identify differences in gene expression that are specific to ivermectin selection and provide a framework to better understand ivermectin resistance and recovery in parasitic helminths. Author Summary Parasitic helminths (worms) infect people and animals throughout the world and are largely controlled with mass administration of anthelmintic drugs. There are a very limited number of anthelmintics available and parasitic helminths can rapidly develop resistance to these drugs. Ivermectin is a widely used anthelmintic in both humans and animals, but resistance is now widespread in the veterinary field. We crossed ivermectin resistant and ivermectin susceptible parasitic helminths and treated them with ivermectin or left them as untreated controls. This provided resistant and susceptible populations with a similar genetic background with which to study differences in gene expression associated with ivermectin resistance. We identified upregulation of a gene with no previous association with drug resistance (HCON_00155390:cky-1) in male and female worms in all resistant populations. This gene is thought to be expressed in the helminth pharynx (mouthpart) and, in mammals, plays a role in controlling nerve function and protecting nerves from damage. This is consistent with the known effects of ivermectin in inhibiting helminth feeding through pharyngeal paralysis and implicates a novel mechanism that allows resistant worms to survive treatment.

populations identified high transcriptomic differentiation between parental isolates, which was significantly reduced in the F2, allowing differences associated specifically with ivermectin resistance to be identified. In all resistant populations, there was constitutive upregulation of a single gene, HCON_00155390:cky-1, a putative pharyngeal-expressed transcription factor, in a narrow locus on chromosome V previously shown to be under ivermectin selection. In addition, we detected sex-specific differences in gene expression between resistant and susceptible populations, including constitutive upregulation of a P-glycoprotein, HCON_00162780:pgp-11, in resistant males only. After ivermectin selection, we identified differential expression of genes with roles in neuronal function and chloride homeostasis, which is consistent with an adaptive response to ivermectin-induced hyperpolarisation of neuromuscular cells. Overall, we show the utility of a genetic cross to identify differences in gene expression that are specific to ivermectin selection and provide a framework to better understand ivermectin resistance and recovery in parasitic helminths.

Author Summary
Parasitic helminths (worms) infect people and animals throughout the world and are largely controlled with mass administration of anthelmintic drugs. There are a very limited number of anthelmintics available and parasitic helminths can rapidly develop resistance to these drugs. Ivermectin is a widely used anthelmintic in both humans and animals, but resistance is now widespread in the veterinary field. We crossed ivermectin resistant and ivermectin susceptible

Introduction
Ivermectin, a 'wonder drug' which won its discoverers the Nobel Prize in Physiology or Medicine, is used primarily to control parasitic disease in humans and animals. Hundreds of millions of people are treated with ivermectin every year to control the filarial nematodes responsible for onchocerciasis (river blindness) and, in combination with albendazole, lymphatic filariasis [1].
Concurrently, ivermectin is one of the top selling animal health products in the world, used extensively in livestock to treat all major gastrointestinal nematodes, in addition to various ectoparasites, and for heartworm prophylaxis in domestic pets. However, ivermectin resistance is now widespread in the veterinary field [2,3] and there are growing concerns of sub-optimal efficacy in humans [4][5][6].
It is hard to overstate the importance of ivermectin in human and animal health, yet there is limited understanding of either the molecular targets of ivermectin in parasitic worms or the mechanisms underlying ivermectin resistance (these may or may not be related). While candidate gene studies have identified single nucleotide polymorphisms (SNPs) or differences in expression of individual genes in resistant and susceptible parasites [7][8][9][10][11], these experiments do not account for differences in the evolutionary histories or genetic backgrounds of the resistant and susceptible populations, so the relevance of these genes to ivermectin resistance remains unclear. Further, exposure to an anthelmintic is likely to have wide-ranging biological effects in parasitic worms, yet the broad scale transcriptional response to drug treatment, and how this translates to the evolution of resistance, is largely unknown.
Haemonchus contortus is a highly pathogenic and economically important gastrointestinal nematode of small ruminants, and ivermectin resistant populations are now common throughout the world [3]. Relative to other parasitic nematodes, H. contortus is a tractable model for anthelmintic resistance research: it resides in the same phylogenetic group as the free-living model nematode Caenorhabditis elegans [12,13] and has a reference quality genome assembly and annotation [14] allowing robust genome-wide analyses [15,16].
To facilitate interrogation of complex traits, such as drug resistance, while controlling for between-population diversity, a genetic cross between genetically divergent populations differing in traits of interest can be generated [15,17]. Progeny of the cross provide a controlled (admixed) population in which to identify genetic variants that co-vary with the trait of interest, for example after drug selection. Applying these principles, we performed a genetic cross between the H. contortus reference genome isolate (MHco3(ISE)) and a multi-drug resistant isolate (MHco18(UGA2004)) and selected adult worms of the F2 generation with ivermectin, benzimidazole or levamisole in vivo [18]. Genomic analyses of the F3 progeny pre-and post-treatment identified discrete loci under selection by each anthelmintic [18] and significantly narrowed a previously identified major quantitative trait locus (QTL) for ivermectin resistance on chromosome V [17].
In this paper, we characterise the broad scale transcriptomic response to ivermectin selection in populations with the same genetic background: male and female adults of the F2 generation of the genetic cross with and without ivermectin treatment. By also measuring differential expression in the two parental isolates, we identify genes specifically associated with ivermectin resistance (i.e. differentially expressed in both the ivermectin-treated genetic cross and the untreated resistant parent) and genes involved in the response to ivermectin treatment (i.e. differentially expressed in the ivermectin-treated genetic cross but not the untreated resistant parent). We highlight a single gene, HCON_00155390:cky-1, in the major QTL for ivermectin resistance, that is constitutively upregulated in all resistant populations, and a P-glycoprotein, HCON_00162780:pgp-11, that is upregulated in resistant males only. Several genes involved in neuronal development and regeneration, neuropeptide signalling and chloride homeostasis were also differentially expressed, providing a framework to better understand the effects of ivermectin on parasitic nematodes.

High transcriptomic differentiation between parental isolates with admixture in the genetic cross
We generated RNA-seq data from male and female worms from both parental isolates (MHco3(ISE); drug susceptible reference isolate, referred to from here on as 'MHco3' and MHco18(UGA2004); multi-drug resistant and referred to from here on as 'MHco18') and from the F2 generation of the genetic cross after ivermectin treatment (F2IVM) and without treatment (F2CTL) ( Figure S1). We also collected samples from the F2 after treatment with benzimidazole (F2BZ) for comparison. An average of 72 million (SD = 27.5 million) 150 bp paired-end reads were generated per sample ( Figure S2A). Between 75 and 82% of MHco18 and MHco3 reads mapped to the reference (MHco3) genome, respectively, with the F2 samples showing intermediate mapping percentages (Figure S2B). This is consistent with the expected levels of polymorphism in each population relative to MHco3, with lower mapping percentages as populations diverge from the reference isolate [18]. Transcriptomic differentiation between isolates, treatment groups, and replicates was assessed using principal component analysis (PCA) and Poisson distance ( Figure 1). Samples clustered well by isolate or treatment group for all male samples ( Figures 1A and 1C); as expected, the parental isolates were highly differentiated (1A, PC1: 38% variance) and the genetic cross F2 samples were an admixture of the two. For the female samples ( Figures 1B and 1D), the parental isolates were again well differentiated, but there was greater variance among all samples (1B, PC1: 53%) with less defined clustering of replicates and treatment groups. Differences between replicates broadly corresponded to the sample processing batch (highlighted by shape of point in 1B and sample name in 1D); this batch effect was included as a secondary factor in the DESeq2 analysis. The female F2CTL group was an outlier in the dataset ( Figure S3), so was excluded from analysis. The specific cause of this variance is unknown: the cluster of genes (n = 129) that differentiated female F2CTL samples and the first replicate of all female samples were mostly uncharacterised; 16 had C. elegans homologues and in this subset there was no enrichment for expression in a particular tissue, or for any GO term. showing Poisson distance between samples. MHco3 = drug-susceptible parent, MHco18 = drug-resistant parent, F2IVM = ivermectin-treated genetic cross F2, F2BZ = benzimidazole-treated F2, F2CTL = untreated control F2. The shape of point in the PCA and the sample number in heatmap relate to the order that the samples were processed (circle = 1st, triangle = 2nd, square = 3rd).

Separating transcriptomic differences associated with drug resistance from other between-isolate variation
To investigate whether ivermectin resistance was associated with differential gene expression, we undertook multiple pairwise comparisons ( Figure S1, Figure   S4). We predicted that the MHco3 and MHco18 isolates would have many constitutive differences in gene expression reflecting genetic diversity resulting from their distinct evolutionary histories, including differences in their sensitivity to the benzimidazoles, levamisole and ivermectin. However, we predicted that strain-specific and drug-specific differences could be separated using the genetic cross: the subset of genes specifically associated with resistance to ivermectin would be differentially expressed in the ivermectin selected genetic cross relative to the unselected genetic cross, but strain-specific differences would not. Consistent with these expectations, the largest number of differentially expressed genes that were unique to a pairwise comparison was between parental isolates (n = 1481 and 1290, males and females respectively, Figure S4). Genes that were only differentially expressed in the parental isolates ( Figure 2, grey points) were essentially ruled out as playing a role in ivermectin resistance mediated by gene expression because they were not also differentially expressed in the ivermectin selected genetic cross ( Figure 2, blue points), with the caveat that our analysis was by choice conservative being based on a stringent cut off (adj P < 0.01) for every pairwise comparison. Using this approach, we identified differentially expressed genes that were specifically associated with ivermectin resistance. comparisons, i.e. likely to represent between-strain differences not associated with ivermectin resistance. If blue they are also differentially expressed in the F2IVM vs MHco3 and F2IVM vs F2CTL comparisons, i.e. likely to be specifically associated with ivermectin resistance. The most highly upregulated gene in resistant isolates is the H. contortus homologue of C. elegans cky-1 (HCON_00155390).

Differentially expressed genes are enriched within the ivermectin QTL on chromosome V after ivermectin selection in the genetic cross
We have previously identified a major effect QTL on chromosome V in response to ivermectin selection in laboratory strains and in field samples of H. contortus [17,18]. We inspected the distribution of differentially expressed genes at this locus and throughout the genome for each pairwise comparison ( Figure S5).
These plots highlighted the large differences in gene expression between parental isolates (S5, A and B) and a reduction in the degree of transcriptomic differentiation in comparisons using the ivermectin treated genetic cross (S5, C, D and E). In comparisons between parental isolates, differentially expressed genes were present throughout the genome with no enrichment for differential expression in the ivermectin resistance locus centred around 37.5 Mb on chromosome V (X 2 = 2.62, P = 0.11). However, the subset of these genes that were also differentially expressed in all comparisons with the genetic cross post ivermectin selection (F2IVM) showed significant enrichment for genes located in the ivermectin resistance locus (X 2 = 37.76, P = 7.23E-10).

Upregulation of a single gene within the ivermectin QTL in all resistant isolates
Seventy-six genes were up-or downregulated in all male and female ivermectin resistant populations (Figure S1B, Figure 2, Table S1). Sixty of the differentially expressed genes had homologues in C. elegans (Table S1) but this gene set showed no tissue, phenotype or GO term enrichment from C. elegans data. The gene with highest upregulation, HCON_00155390, a homologue of C. elegans cky-1, lies within the major chromosome V locus under ivermectin selection in H. contortus ( Figure 3). As shown in Table S1, HCON_00155390:cky-1 is expressed between 62 and 132-fold higher (log2 fold changes from 5.95 to 7.04, adj P ≤ with RT-qPCR [18]. Constitutive upregulation in MHco18 females was also detected relative to MHco3 females, but the fold-change in expression could not be calculated due to a lack of measurable expression in MHco3 females (data not shown). In C. elegans, cky-1 encodes a transcription factor expressed in the embryo and in non-neuronal pharyngeal cells [19] and the only published phenotype is early onset of polyglutamine-YFP aggregation after RNAi [20]. The human orthologue, NPAS4 or NXF, is a basic helix-loop-helix transcription factor that regulates gamma-aminobutyric acid (GABA) releasing inhibitory synapses [21] and can be induced in response to neurodegeneration or excitation to confer protection to neuronal cells [22]. Notably, no putative ivermectin target or resistance genes from the literature [7,8,11,23,24] were identified as differentially expressed in all resistant populations; normalised read counts for candidate genes are shown in Figures   S6 (males) and S7 (females). Genes that are differentially expressed in the parental isolates only are unlikely to be related to ivermectin resistance and none of these candidate genes lie within genomic loci under ivermectin selection [17,18].

Male and female worms show shared and sex-specific differences in gene expression associated with ivermectin resistance
We were concerned that the higher between-replicate variability of the female data could limit the detection of truly differentially expressed genes in the male analysis. It is also possible that ivermectin affects male and female worms differently [25,26]. For these reasons, we also analysed results from the male and female datasets separately (Figure S1B, Figure 4). In the male worms, there were an additional 67 upregulated and 107 downregulated genes that were not differentially expressed in all pairwise comparisons of ivermectin resistant female worms (Table S2). There was low but significant upregulation of two putative multi-drug resistance proteins associated with ivermectin resistance in the literature: HCON_00162780 and HCON_00175410, the homologues of C. elegans pgp-11 and mrp-4, respectively. PGP-11 is a P-glycoprotein known to modulate ivermectin sensitivity in other nematode species [27], and the pgp-11 locus has previously been associated with ivermectin resistance in H. contortus [16,17,24]. In this study, the expression of HCON_00162780:pgp-11 was markedly higher in males relative to females ( Figure 5A and 5B), and upregulation of pgp-11 in MHco18 versus MHco3 males was confirmed with RT-qPCR ( Figure 5C). In C. elegans, mrp-4 is involved in lipid transport and lipid storage, but has also been shown to have higher constitutive expression in an ivermectin resistant strain of C. elegans relative to wild-type [28]. One putative ligand gated ion channel, HCON_00137455, a homologue of lgc-25, was downregulated, as was HCON_00078660, a homologue of the sodium leak ion channel nca-2. In C. elegans, nca-2 is expressed in acetylcholine and GABA motor neurons [29], and is involved in the propagation of neuronal activity from cell bodies to synapses [30] and in synaptic vesicle recycling [29]. There was upregulation of one chloride channel, HCON_00108990, the homologue of best- 19. In C. elegans, this is a bestrophin chloride channel expressed in the NSM, a pharyngeal neurosecretory motor neuron.  Of the differentially expressed genes in ivermectin resistant males, 125 had C. elegans homologues (Table S2). These homologues were enriched for expression in the ALA neuron (E = 0.24, O = 4, Q = 0.0011) ( Figure S8), a mechanosensory interneuron that projects into the nerve ring then joins the lateral cords, extending as far as the tail, adjacent to the excretory canals [31]. The ALA neuron inhibits locomotion and pharyngeal pumping during lethargus and 'stress induced sleep', a characteristic behaviour in C. elegans during recovery from exposure to damaging conditions [32,33]. There was GO enrichment for a single Notably, there was also phenotype enrichment for 13 terms associated with movement ( Figure S8).
For female worms, the lack of an F2CTL sample meant the equivalent analysis (excluding differentially expressed genes in all resistant male comparisons) would be very likely to be confounded by differences between the female parental isolates that were not associated with ivermectin resistance ( Figure   4B). To mitigate this, an alternative analysis including previously published RNAseq data for females from two additional ivermectin resistant populations (MHco4(WRS) and MHco10(CAVR) versus MHco3 [34]) was undertaken ( Figure   S1B). This identified 85 upregulated and 97 downregulated genes in all ivermectin resistant female samples that were not differentially expressed in ivermectin resistant male samples (Table S3). Downregulation of HCON_00076370, a homologue of the C. elegans sodium:potassium:chloride symporter nkcc-1, was identified in this analysis. In C. elegans, nkcc-1 is a chloride transporter and functions as a developmental switch that changes the response of GABAA receptors in body wall muscles from depolarising in the early L1 to hyperpolarising from the mid-L1 stage onwards by controlling the cellular chloride gradient [35]. In female worms only, there was also downregulation of one putative ligand gated ion channel, HCON_00097440, a homologue of acr-24.
Of the differentially expressed genes in ivermectin resistant females, 158 had C. elegans homologues. These showed no tissue or GO term enrichment but there was enrichment for a single phenotype 'uptake by intestinal cell defective' (E = 0.59, O = 4, Q = 0.075).

Transcriptomic response to ivermectin exposure implicates neuronal regeneration, neuropeptide signalling and chloride homeostasis
Having identified genes with constitutive differential expression in ivermectin resistant worms, we were also interested in genes that were only up or downregulated after exposure to ivermectin treatment. This subset of genes would be expected to be differentially expressed in the F2IVM versus F2CTL comparison and the F2IVM versus MHco3 comparisons, but not the MHco18 versus MHco3 comparisons, because the MHco18 isolate had not been exposed to ivermectin. As ivermectin has a short half-life in sheep (plasma t1/2 = 2.7 days [36]) and the F2IVM worms were isolated 21 days after ivermectin treatment, this gene set will not reflect acute ivermectin exposure but rather the longerterm influence of drug treatment.
Of 120 genes that were differentially expressed in response to ivermectin treatment, 95 had C. elegans homologues (Table S4). The most highly upregulated gene with a C. elegans homologue was HCON_00068595 which was 5 to 32-fold more highly expressed in ivermectin treated populations (2.31 -4.99 log2 fold change, adj P ≤ 2.52E-12 (Table S4)). This gene is a homologue of emb-9, a collagen type IV alpha 5 chain, involved in many processes such as pharyngeal pumping and egg laying, but with a critical role in positive regulation of neuron regeneration following axon damage [37]. emb-9 is most strongly expressed in the basement membrane of the C. elegans pharynx [38]. There was no differential expression of any ligand gated ion channels after ivermectin exposure, but HCON_00068600, the homologue of a neuronal CLC-2 chloride channel, clh-3, was upregulated between 2 and 12-fold (1.28 -3.62 log2 fold change, P ≤ 2.95E-10 (Table S4)

Differential expression of genes associated with neuronal plasticity and chloride homeostasis are specific to ivermectin selection
To investigate the specificity of the transcriptomic differences we identified in the ivermectin-selected genetic cross and to explore the possibility of a common response to drug exposure (for example, a stress response or multi-drug resistance pathway), we also undertook RNA-seq analysis of the same F2 population after benzimidazole selection. There were fewer differentially expressed genes in the benzimidazole resistant populations relative to the equivalent comparisons for ivermectin (59 associated with benzimidazole resistance (Table S5) versus 76 for ivermectin, and 55 in response to benzimidazole treatment (Table S6) versus 120 for ivermectin). The benzimidazole resistant populations showed no differential expression of genes known to confer resistance through non-synonymous mutations in the target protein, -tubulin (i.e., -tubulin isotype 1 (HCON_00005260) or -tubulin isotype 2 (HCON_00043670)) or genes associated with drug metabolism (UDPglycotransferases or cytochrome P450 enzymes).
Notably, both the ivermectin and benzimidazole-selected populations appeared to be transcriptionally very similar (Figure 1), despite different regions of the genome being under selection by each drug, with no evidence of a shared genetic component to resistance [18]. In the benzimidazole resistant and response to treatment groups, 35 genes in each (59% and 65% of differentially expressed genes, respectively) were shared with the equivalent ivermectin selected However, in contrast with the ivermectin resistant populations, there was no enrichment for differential expression at the chromosome V locus (Figure 3, bottom panel) and no differential expression of HCON_00155390:cky-1 or any of the genes associated with neuronal regeneration and chloride homeostasis. This suggests the differential expression of these genes, identified using the same genetic cross, are unique to ivermectin resistance and recovery.

Use of a genetic cross minimises the impact of between-isolate polymorphism on differential expression analysis
Previous work highlighted the potential impact of high levels of polymorphism on between-isolate RNA-seq analysis in H. contortus [34], showing that polymorphic isolates (relative to the reference MHco3 isolate) had significantly fewer reads mapping to the genome, potentially biasing differential expression analysis. In this study, the difference in 'SNP rate' (the number of SNPs in a gene model; see Methods) between parental isolates had a small but significant negative correlation with gene expression (R 2 = 0.074 [P < 0.001] and R 2 = 0.099 [P < 0.001] for males and females, respectively). This correlation suggests that highly polymorphic genes were more likely to be identified as downregulated than upregulated in the parental isolates ( Figure S9). This effect was more striking for the MHco18 isolate (upper left quadrant) but could be seen to a lesser extent for polymorphic genes in the MHco3 isolate (lower right quadrant).

Discussion
A key challenge of genomic analyses of helminths is accounting for the high genetic diversity between populations [40]. Comparisons between resistant and susceptible populations are further complicated by inherent differences in traits other than resistance [16]. In the case of transcriptomic analyses, the same factors confound detection of differential gene expression specifically associated with resistance, and commonly used measurements of gene expression can be biased by between-population polymorphism [34]. Further, for transcriptomic analyses, differential gene expression after drug selection may reflect constitutive differences in gene expression between resistant and susceptible worms and/or differences that are induced by drug exposure. The latter may be specific to a particular drug or may represent a more generic stress response or drug metabolism pathway. Separating these different mechanisms in betweenpopulation comparisons is complex and usually not feasible.
In this study, we used a genetic cross of a multi-drug resistant isolate of H. contortus with a susceptible (reference) isolate, providing a population composed of an admixture of parental genotypes, to control for differences in the genetic background of each parent [18] and to mitigate the impact of polymorphism on downstream analyses. Comparisons of resistant adults with and without ivermectin selection allowed us to separate differential expression associated with ivermectin resistance from that associated with the response to drug exposure. Further, a comparison of genes with differential expression after selection with ivermectin or benzimidazole allowed us to identify drug-specific and shared responses.
The success of the genetic crossing approach in reducing the impact of background genetic diversity between the parental isolates was demonstrated by the reduction in the number of differentially expressed genes identified in the genetic cross with and without ivermectin treatment, relative to the number of differentially expressed genes identified between parents. In the genetic cross, there was also a reduction in the distribution of differentially expressed genes throughout the genome, with significant enrichment for differentially expressed genes within the ivermectin resistance QTL on chromosome V, which was not detectable in the parental comparison due to high transcriptomic diversity throughout the genome.
A previous study highlighted the potential confounding effect of genetic diversity on RNA-seq analysis of different H. contortus isolates due to poorer mapping of divergent reads to the reference genome [34]. In this study, alignment of longer reads (150 bp vs 100 bp) to a more contiguous genome assembly [14] with an improved aligner (HISAT2 [41]) resulted in more reads from the divergent MHco18 isolate aligned to the reference genome (75-77%), than from the reference isolate in the previous study (69%) [34]. However, the small but significant negative correlation between level of polymorphism and differential expression in the parental isolates suggests that mapping biases could still impact the accuracy of differential expression measurement for polymorphic genes in between-isolate comparisons. Further, the significant enrichment for an 'anthelmintic response variant' phenotype in the most highly polymorphic subset of genes highlights a challenge of studying complex traits in rapidly evolving helminths: the polymorphic genes may be those of most interest, but their analysis can be confounded when relying on a single reference genome for read alignment and/or primer design, or when comparing sequences from different populations to infer a relationship with the trait of interest. Improved tools for mapping RNA-seq reads to a reference genome, such as HISAT-genotype [42], which incorporates known variants to optimise the mapping of polymorphic regions of the human genome is one promising approach, although such an approach is not yet available for non-model organisms. However, in this study, the use of the genetic cross resulting in an admixture of parental genotypes largely corrected for artefactual differential expression due to mapping bias in male samples, although some bias remained for females due to the lack of a F2 control.
Notably, after controlling for differences in the genetic background of our resistant and susceptible isolates with the genetic cross, we found no evidence of differential expression of candidate ivermectin resistance genes from the literature [7,8,11,23,24], with the exception of HCON_00162780:pgp-11, which was upregulated in resistant males only. While this does not rule out a role for differential expression of these genes in the early stages of resistance development, it suggests they do not confer high level resistance; this is consistent with our genomic analyses of ivermectin resistance, where we find no candidate resistance genes in the major locus under ivermectin selection [17,18]. Despite finding no evidence of differential expression of any glutamate gated chloride channel receptor subunits, the predicted major targets of ivermectin [43][44][45], we did identify differential expression of multiple neuronal genes predicted to modulate cellular chloride levels (clh-3, best-19, nkcc-1), which could potentially mitigate the effects of chloride influx when ivermectin binds to its receptor.
The gene with the most robust association with ivermectin resistance was the homologue of C. elegans cky-1, a gene with no previous association with drug resistance. HCON_00155390:cky-1 lies at the centre of the major locus of ivermectin resistance on chromosome V [18] and was the most highly upregulated gene overall in both males and females. In C. elegans, cky-1 is a transcription factor encoding basic helix-loop-helix (bHLH) and Per-Arnt-Sim (PAS) domains; these domains are also present in the H. contortus homologue.
Very little is known about the function of cky-1 in nematodes, but the mammalian orthologue (NPAS4/NXF) regulates GABA releasing inhibitory synapses [21] and is induced by neurodegeneration or excitation to confer protection to neuronal cells [22]. There are many hundred predicted targets of NPAS4 in the mouse, based on differential expression measured by microarray after NPAS4-RNAi [21], including various ion channels and synaptic proteins, but many (more than a third) of the target genes are uncharacterised. In C. elegans, only four cky-1-interacting genes are listed in WormBase: aha-1, pmk-1, vab-3 and ztf-2. Homologues of all four genes are present in the H. contortus genome but none are differentially expressed in all ivermectin resistant adults. It is possible that some of the other differentially expressed genes identified in this study are downstream targets of HCON_00155390:cky-1 and this warrants further investigation. In C. elegans, cky-1 expression is highest in the embryo and is largely restricted to pharyngeal cells (arcade cells, pharyngeal muscle cells, hypodermis, and pharyngeal:intestinal valve) with little or no expression detected elsewhere by single cell RNA-seq [19]. It is unknown if this pharyngeal expression pattern is the same in either ivermectin susceptible or resistant H. contortus but studies are underway. The major routes of ivermectin uptake in adult H. contortus are unclear, but it is likely that the pharynx plays a key role.
Despite surviving treatment, resistant nematodes can show phenotypic effects after ivermectin exposure, for example, reduced pumping in the triple resistant C. elegans DA1316 [46] and suppression of egg production in gastrointestinal nematodes [47], particularly at the early stages of resistance development, which suggests they are still physiologically affected to some degree by the drug. These phenotypic effects could relate to the many potential targets of ivermectin [48] and/or the mode of resistance. We identified significant enrichment for downregulation of multiple inhibitory neuropeptides after drug exposure (and with constitutive differential expression in resistant males). These signalling peptides regulate various aspects of behaviour, including locomotion, pharyngeal pumping, lipid storage and egg laying, which may relate to phenotypes observed in ivermectin treated worms. Transcriptomic data also highlighted genes associated with neuronal plasticity, neuronal development and regeneration, yet many of these genes are pleiotropic, and differential expression may equally relate to their roles in overcoming the phenotypic effects of ivermectin on feeding, movement and reproduction. For example, in C. elegans, emb-9 positively regulates neuronal regeneration following damage yet is also essential for pharyngeal pumping and egg laying [37].
Our analysis of the same genetic cross after benzimidazole selection supports our expectation that the genes we have identified with putative roles in ivermectin resistance are largely unique to this drug. However, the degree of overlap in the two groups was higher than expected given the separate modes of action and distinct genomic loci under selection by each drug [18]. One technical explanation could be the impact of the male F2CTL sample: if, as seen in the female data, this sample was an outlier, it could artificially inflate the number of shared differentially expressed genes in the F2BZ and F2IVM samples. Despite passing all our QC analysis, there was some evidence to support this in the higher number of genes that were differentially expressed in the male F2IVM versus F2CTL comparison relative to the male F2IVM versus MHco3 comparison, which is hard to explain biologically. However, even if this is the case, based on the PCA plots, the benzimidazole-selected and ivermectin-selected populations do appear to be transcriptomically similar. Another explanation would be stochastic inheritance of the same parental alleles at regions of the genome that are not necessarily under selection by both drugs; this would reflect the divergence of the two parental isolates and is increasingly likely with a small sample size (in this case, 20 worms per sample), in a polyandrous species [49], and with a limited number of recombination events (F2 generation of the cross).
Alternatively, the results could indicate a degree of co-selection of mechanisms promoting tolerance of multiple anthelmintics in a population routinely treated with all drug classes. However, other than pgp-11, which may represent a multidrug resistant protein, there were no known resistance-associated genes or pathways shared between the ivermectin and benzimidazole selected groups.
One particularly striking finding was the apparent lack of differential expression of genes involved in xenobiotic metabolism for either drug class; this is consistent with the expectation that ivermectin is not metabolised in nematodes [46], but is more surprising for the benzimidazoles where biotransformation has been detected in H. contortus and C. elegans [50][51][52]. This may reflect the extended time between treatment and isolation of worms in this study, which would miss an acute response to circulating drug. However, our results are consistent with a recent study investigating the acute transcriptional response of adult MHco18 H. contortus exposed to different benzimidazoles in vitro, where, in contrast to the equivalent experiment in C. elegans, no genes with roles in xenobiotic metabolism were differentially expressed [53].
In conclusion, the use of a genetic cross has allowed the separation of transcriptomic differences associated with ivermectin resistance and survival from a background of high transcriptomic differentiation in resistant and susceptible isolates of H. contortus. We identify constitutive upregulation of HCON_00155390:cky-1, a gene in the centre of the major genomic locus under ivermectin selection in global populations of H. contortus. We also identify a small number of differentially expressed genes outwith the genomic locus, which may be downstream targets of cky-1 or may function in shared or parallel pathways to promote ivermectin resistance.

Generation of parasite material
Fifteen donor sheep were treated with 2.5 mg/kg bodyweight monepantel (Zolvix Oral Solution for Sheep, Elanco AH) 14 days prior to infection to ensure they were parasite free. Each donor was orally infected with 5000 L3 of a single H. contortus isolate; three donors were infected with the drug susceptible MHco3(ISE) isolate [54], three donors were infected with the triple resistant MHco18(UGA2004) isolate (benzimidazole, levamisole and ivermectin resistant) [7] and nine donors were infected with the F2 generation of a genetic cross between MHco3 females and MHco18 males [15].

Gene enrichment analysis
Putative C. elegans homologues were identified with BLASTP (C. elegans WS264 proteins, e-value < 0.001) and tissue, phenotype and gene ontology enrichment analysis was performed with WormBase Enrichment Analysis [68], which uses a hypergeometric test and Benjamini-Hochberg correction of false discovery rate (Q < 0.1 for significance).

Real-time quantitative PCR (RT-qPCR)
Total RNA was extracted from 20 male or female worms from each of three different donor sheep per isolate (MHco3 and MHco18). Total RNA (1 μg) was used for oligo(dT) primed cDNA synthesis (SuperScript® III First-Strand Synthesis System, ThermoFisher, 18080051), with no-reverse transcriptase controls included for each sample. cDNA was diluted 1:10 for RT-qPCR (other than for HCON_00155390:cky-1 in female samples, where neat cDNA was used due to the low abundance of this transcript) and 1 μl added to each reaction. RT-qPCR was undertaken following the Brilliant III Ultra Fast SYBR QPCR Master Mix protocol (Agilent Technologies, 600882) and results were analysed with MxPro v4.10. Gene expression was normalised to β-actin (HCON_00135080).
Primer sequences are listed in Table S7.  S1. Study design. A. The Haemonchus contortus MHco3 isolate is fully drug susceptible and the MHco18 isolate is multi-drug resistant; these are the two parental populations used to generate the genetic cross, of which the parent and F2 adult populations are used in this study. One donor sheep was infected per population. Nomenclature, drug treatment, expected resistance phenotype(s) and adult worm samples used in our analyses are described. B. Pairwise comparisons of differential gene expression included in each analysis. For each analysis, genes were filtered for significant differential expressed (adjusted P<0.01) in all pairwise comparisons marked ✓ and for non-significant differential expression (adjusted P>0.01) in all pairwise comparisons marked ✕.

Acknowledgements
Pairwise comparisons marked -were not included in the analysis. Genes with significant differential expression were filtered to retain only those where the direction of log fold change (>0 or <0) was the same in all pairwise comparisons i.e. for inclusion, a gene could not be upregulated in a resistant population in one pairwise comparison and downregulated in a resistant population in another.