(C) PLOS One
This story was originally published by PLOS One and is unaltered.
 . . . . . . . . . .



Transcriptomic analyses implicate neuronal plasticity and chloride homeostasis in ivermectin resistance and response to treatment in a parasitic nematode [1]

['Roz Laing', 'Institute Of Biodiversity Animal Health', 'Comparative Medicine', 'College Of Medical', 'Veterinary', 'Life Sciences', 'University Of Glasgow', 'Glasgow', 'United Kingdom', 'Stephen R. Doyle']

Date: 2022-08 

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 response to treatment in parasitic helminths.

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 may implicate a novel mechanism that allows resistant worms to survive treatment.

Funding: This work was funded by a Biotechnology and Biological Sciences Research Council (BBSRC) strategic Lola [BB/M003949] (RL, DJB, NS, JAC, CB, ED), Scottish Government’s Rural and Environment Science and Analytical Services (RESAS) Division (DJB) and by the Wellcome Trust [206194]. RL is supported by a Wellcome Clinical Research Career Development Fellowship [216614/Z/19/Z], SRD is supported by a UKRI Future Leaders Fellowship [MR/T020733/1]. The funders had no role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript. For the purpose of Open Access, the author has applied a CC BY public copyright licence to any Author Accepted Manuscript version arising from this submission.

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.

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 treated the F2 generation of adult worms with standard doses of ivermectin, benzimidazole or levamisole in vivo [ 18 ]. Genomic analyses of the F3 progeny collected pre- and post-treatment of the F2 adults that survived 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 ].

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 – 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.

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 – 6 ].

The positive impact of the cross in the male analyses was further supported by a reduction in the proportion of highly polymorphic genes (defined as having a SNP rate of greater than 2% in MHco18) in the downregulated genes identified in the F2IVM versus F2CTL analysis compared to the MHco18 versus MHco3 analysis. Notably, 21% (469/2196) of genes identified as downregulated in MHco18 males were highly polymorphic, relative to 2.8% (58/2069) of upregulated genes, whereas in the genetic cross 6.2% (60/973) of genes identified as downregulated in F2IVM males were highly polymorphic, relative to 5.8% (71/1219) of upregulated genes. Finally, the relevance of controlling for polymorphism was highlighted by the subset of genes that were identified as the most highly polymorphic in the parental isolates: C. elegans homologues of the 1684 genes with a greater than 2% difference in polymorphism rate in the MHco18 isolate relative to the MHco3, were enriched for the phenotype ‘anthelmintic response variant’ (E = 3.2, O = 10, Q = 0.058). This was associated with highly polymorphic homologues of avr-15, acr-16, lev-1, osm-1, che-1, che-3, che-11, che-14, daf-10 and dyf-3. Although not identified in the enrichment analysis, a homologue of ben-1, HCON_00043670 (β-tubulin isotype 2), was also highly polymorphic. Overall, these findings highlight the need to control for polymorphism in transcriptomic analyses of genetically divergent populations and show that genes associated with drug response are among the most polymorphic genes in a resistant population.

Previous work highlighted the potential impact of high levels of polymorphism on between-isolate RNA-seq analysis in H. contortus [ 30 ], 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 Materials and 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 ( S10 Fig ). 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). Admixture of the parental genotypes in the genetic cross largely mitigated this bias in male analyses (R 2 = 0.0026 [P = 0.322]) ( S10 Fig , dark blue points), although the lack of a female F2CTL limited our ability to control for polymorphism in the female-only analyses (R 2 = 0.074 [P < 0.001]).

However, in contrast with the ivermectin resistant populations, there was no enrichment for differential expression at the chromosome V locus ( Fig 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 response to treatment.

Notably, both the ivermectin and benzimidazole-selected populations appeared to be transcriptionally very similar ( Fig 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 populations. In the benzimidazole resistant populations, as for the ivermectin resistant populations, there was downregulation of multiple homologues of C. elegans vap-1 (five downregulated constitutively and 11 in response to benzimidazole treatment) and four of the same neuropeptides were downregulated in response to benzimidazole and ivermectin treatment (homologues of nlp-13, flp-1, flp-14 and flp-22). In male worms only, upregulation of one putative multi-drug resistant gene, HCON_00162780:pgp-11, was identified in all pairwise comparisons for both ivermectin and benzimidazole resistant populations.

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 ( S5 Table ) versus 76 for ivermectin, and 55 in response to benzimidazole treatment ( S6 Table ) 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 (UDP-glycotransferases or cytochrome P450 enzymes).

For the C. elegans homologues, the most significant GO term was ‘neuropeptide signalling pathway’ (E = 0.71, O = 8, Q = 6e-06) ( S9 Fig ). The eight genes associated with this term encoded five neuropeptide-like proteins (nlp-7, nlp-11, nlp-13, nlp-49 and nlp-50) and three FMRF-like peptides (flp-1, flp-14 and flp-22); all were downregulated in F2IVM samples. In common with the constitutive datasets, multiple homologues of C. elegans vap-1 (15 in total) were downregulated.

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 ( S4 Table )) in ivermectin treated populations. CLC-2 is activated by membrane hyperpolarisation and is thought to regulate neuronal excitability by providing an efflux pathway for chloride. In C. elegans, clh-3 is only expressed in the (spontaneously active) hermaphrodite-specific motor neuron (HSN), where clh-3 inhibits excitability by mediating chloride influx, resulting in suppression of egg laying [ 35 ]. In H. contortus, HCON_00068600:clh-3 is also constitutively expressed in both male parental samples, suggesting expression is not limited to the HSN in this species. Interestingly, HCON_00068595:emb-9 and HCON_00068600:clh-3 are neighbouring genes on chromosome III, although encoded on opposite strands. In C. elegans, emb-9 and clh-3 are encoded on different chromosomes (chromosome III and II, respectively) and have no known interaction.

Of 120 genes that were differentially expressed in response to ivermectin treatment, 95 had C. elegans homologues ( S4 Table ); none of these genes were within the major QTL. 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 ( S4 Table )). 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 [ 33 ]. emb-9 is most strongly expressed in the basement membrane of the C. elegans pharynx [ 34 ].

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 t 1/2 = 2.7 days [ 32 ]) and the F2IVM worms were isolated 21 days after ivermectin treatment, this gene set will not reflect acute ivermectin exposure but rather the longer-term influence of drug treatment.

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 ( Fig 4B ). To mitigate this, an alternative analysis including previously published RNA-seq data for females from two additional ivermectin resistant populations (MHco4(WRS) and MHco10(CAVR) versus MHco3 [ 30 ]) was undertaken ( S1B Fig ). This identified 85 upregulated and 97 downregulated genes in all ivermectin resistant female samples that were not differentially expressed in ivermectin resistant male samples ( S3 Table ). 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 GABA A 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 [ 31 ]. 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).

Of the differentially expressed genes in ivermectin resistant males, 125 had C. elegans homologues ( S2 Table ). These homologues were enriched for expression in the ALA neuron (E = 0.24, O = 4, Q = 0.0011) ( S9 Fig ), 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 [ 27 ]. 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 [ 28 , 29 ]. There was GO enrichment for a single term ‘neuropeptide signalling pathway’ (E = 0.92, O = 6, Q = 0.013), associated with nlp-9, nlp-14, flp-6, flp-7, flp-16 and flp-24, all of which were downregulated. Notably, there was also phenotype enrichment for 13 terms associated with movement ( S9 Fig ).

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 when results from both datasets were combined. It is also possible that ivermectin affects male and female worms differently [ 21 , 22 ]. For these reasons, we also analysed results from the male and female datasets separately (Figs 4 and S1B ). None of the genes identified in these sex-specific analyses were within the major QTL. 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 ( S2 Table ). 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 [ 23 ], and the pgp-11 locus has previously been associated with ivermectin resistance in H. contortus [ 16 , 17 , 20 ]. In this study, the expression of HCON_00162780:pgp-11 was markedly higher in males relative to females ( Fig 5A and 5B ), and upregulation in MHco18 versus MHco3 males was confirmed with RT-qPCR ( Fig 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 [ 24 ]. 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 [ 25 ], and is involved in the propagation of neuronal activity from cell bodies to synapses [ 26 ] and in synaptic vesicle recycling [ 25 ]. 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.

Notably, no putative ivermectin target or resistance genes from the literature [ 7 , 8 , 11 , 20 ] were identified as differentially expressed in all resistant populations; normalised read counts for candidate genes are shown in S7 (males) and S8 (females) Figs. 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 ].

One of the most highly downregulated genes also lies within the major locus under ivermectin selection: HCON_00155240, a homologue of C. elegans F59B1.10. In C. elegans, this gene encodes a checkpoint (CHK) kinase-like protein with an uncharacterised oxidoreductase Dhs-27 domain, expressed in GABAnergic neurons, dopaminergic neurons and muscle cells. One caveat to our measurement of expression for this gene is its high level of polymorphism in the MHco18 isolate, which could give the appearance of downregulation due to a greater proportion of unmapped reads relative to MHco3 (discussed in more detail below). Outside of the ivermectin QTL, downregulated genes included seven homologues of C. elegans vap-1, a cysteine rich secretory protein expressed in the amphid sheath cell.

Chromosome V karyoplots showing genomic position of genes with significant upregulation (yellow) or downregulation (blue). Each point represents a differentially expressed gene (adj P < 0.01) and point size corresponds to significance. Panels A-D compare ivermectin-resistant versus ivermectin-susceptible populations (A. MHco18 vs MHco3 males, B. MHco18 vs MHco3 females, C. F2IVM vs MHco3 males, D. F2IVM vs MHco3 females). Panel E compares an ivermectin-resistant versus a mixed (resistant and susceptible) population (F2IVM vs F2CTL males). Panel F shows genetic differentiation (F ST ) between the F3 generation of the genetic cross pre- and post- ivermectin selection to highlight the major locus of ivermectin selection [ 18 ]. Panel G shows relative absence of differential expression at the ivermectin QTL in a comparison with a benzimidazole selected population (F2BZ vs F2CTL males). Upregulation of the H. contortus homologue of C. elegans cky-1 (HCON_00155390) is identified in ivermectin resistant populations only.

Seventy-six genes were up- or downregulated in all male and female ivermectin resistant populations (Figs 2 and S1B and S1 Table ). Sixty of the differentially expressed genes had homologues in C. elegans ( S1 Table ) 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 ( Fig 3 ). As shown in S1 Table , HCON_00155390:cky-1 is expressed between 62 and 132-fold higher (log2 fold changes from 5.95 to 7.04, adj P ≤ 0.00024) in males and females of all resistant populations relative to the susceptible MHco3 isolate and 2.7-fold higher (1.44 log2 fold change, adj P = 0.0019) in males of the genetic cross after ivermectin selection (F2IVM = resistant genotypes only) relative to the unselected population (F2CTL = mix of resistant and susceptible genotypes). In the latter comparison the control population is an admixture of resistant and susceptible genotypes so the relative difference is expected to be lower than in the former comparisons with the fully susceptible MHco3 control. Constitutive upregulation of HCON_00155390:cky-1 was confirmed in MHco18 males relative to MHco3 males 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. In C. elegans, cky-1 encodes a transcription factor expressed in the embryo and in non-neuronal pharyngeal cells [ 19 ].

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 ( S6 Fig ). These plots highlighted the large differences in gene expression between parental isolates (S6, A and B) and a reduction in the degree of transcriptomic differentiation in comparisons using the ivermectin treated genetic cross (S6, 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).

Scatter plot showing differentially expressed genes in comparisons of resistant and susceptible populations (adj P < 0.01), with log2 fold change of males shown on the x-axis and females on the y-axis. Grey points represent genes that are differentially expressed only in parental (MHco18 vs MHco3) 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). Equivalent scatterplots for all pairwise comparisons are shown in S5 Fig .

To investigate whether ivermectin resistance was associated with differential gene expression, we undertook multiple pairwise comparisons ( S1 and S4 Figs). 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, S4 Fig ). Genes that were only differentially expressed in the parental isolates ( Fig 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 ( Fig 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.

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) ( S1 Fig ). 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 ( S2A Fig ). Between 75 and 82% of MHco18 and MHco3 reads mapped to the reference (MHco3) genome, respectively, with the F2 samples showing intermediate mapping percentages ( S2B Fig ). 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 ( Fig 1 ). Samples clustered well by isolate or treatment group for all male samples ( Fig 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 ( Fig 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 datapoint symbol 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 ( S3 Fig ), 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.

Discussion

A key challenge of genomic analyses of helminths is accounting for the high genetic diversity between populations [36]. 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 [30]. 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 between-population 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 [30]. In this study, alignment of longer reads (150 bp vs 100 bp) to a more contiguous genome assembly [14] with an improved aligner (HISAT2 [37]) 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%) [30]. 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 [38], 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,20], 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 [39–41], 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 [42] and is induced by neurodegeneration or excitation to confer protection to neuronal cells [43]. There are hundreds of predicted targets of NPAS4 in the mouse, based on differential expression measured by microarray after NPAS4-RNAi [42], 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. The major routes of ivermectin uptake in adult H. contortus are unclear, but it is likely that the pharynx plays a key role.

If overexpression of cky-1 confers resistance, then we hypothesise that knock down of cky-1 in a resistant isolate should restore sensitivity to ivermectin. However, to date, we have been unable to knock down HCON_00155390:cky-1 in Haemonchus by RNAi, despite significant effort. Although several previous studies have reported successful gene knock down in Haemonchus L1, L2 and L3 by soaking and/or feeding with dsRNA or siRNAs [44–47], these methods do not work reliably for all genes tested in this species [48,49] unlike in C. elegans. This may relate to intrinsic differences in mechanisms of gene silencing between free-living and parasitic nematodes, possibly due to inefficient uptake and/or amplification of siRNAs in the latter. In addition, phenotypes are challenging to monitor in parasitic nematodes. We have, however, demonstrated increased ivermectin sensitivity in a C. elegans cky-1 loss of function mutant and by RNAi in N2 worms using a progeny survival assay [18]. These experiments show that altered expression of cky-1 can modulate survival of nematodes exposed to ivermectin. While the data from C. elegans are compelling and highlight the strength of this model system to study less tractable parasitic species [50], these data assume that function is conserved between the orthologues of cky-1 in both species. Future work to test for functional conservation by transgenic rescue of lethality in C. elegans cky-1 loss of function mutants and to establish cky-1 RNAi in H. contortus are required to understand if cky-1 plays a role in ivermectin resistance in parasitic worms.

Despite surviving treatment, resistant nematodes can show phenotypic effects after ivermectin exposure, for example, reduced pumping in the triple resistant C. elegans DA1316 [51] and suppression of egg production in gastrointestinal nematodes [52], 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 [53] 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 [33].

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 [54], 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 multi-drug 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 [51], but is more surprising for the benzimidazoles where biotransformation has been detected in H. contortus and C. elegans [55–57]. 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 [58].

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.

[END]
---
[1] Url: https://journals.plos.org/plospathogens/article?id=10.1371/journal.ppat.1010545

Published and (C) by PLOS One
Content appears here under this condition or license: Creative Commons - Attribution BY 4.0.

via Magical.Fish Gopher News Feeds:
gopher://magical.fish/1/feeds/news/plosone/