Comparative Analysis on Transcriptomics of Ivermectin Resistant and Sensitive Strains of Haemonchus Contortus


 Background: Ivermectin (IVM) is one of the most important and widely used anthelmintics in veterinary medicine. However, its efficacy is increasingly compromised by widespread resistance, and the exact mechanism of IVM resistance remains unclear for most parasitic nematodes including Haemonchus contortus, a blood-sucking parasitic nematode of small ruminants.Methods: In this study, we isolated and assessed an IVM resistant strain from Zhaosu, Xinjiang, China. Subsequently, the comparative analyses on transcriptomics of IVM susceptible and resistant H. contortus adult worms were carried out using RNA sequencing and bioinformatics.Results: In total, 543 and 359 differentially expressed genes (DEGs) were identified in male and female adult worms of the resistant strain compared with the susceptible strain, respectively. The DEGs encode molecules involved in receptor activities, transport, detoxification, lipid metabolism and cuticle collagen formation. In addition, Gene Ontology (GO) analysis revealed that transcriptional changes were dominant in genes associated with ligand-gated channel activity, oxidation-reduction process, lipid metabolic process, and structural constituent of cuticle. The results support previous proposal that the IVM resistant mechanism of H. contortus involved in both neuromuscular and non-neuromuscular pathways. Finally, the quantitative RT-PCR results confirmed that the transcriptional profiles of selected DEGs (male: 8 genes, female: 10 genes) were consistent with those obtained by the RNA-Seq.Conclusions: The findings from this work provided valuable information for further studies on the IVM resistance in H. contortus.


Background
Haemonchus contortus is the most prevalent and pathogenic gastrointestinal nematode (GIN) of small ruminants, causing signi cant economic impacts on livestock production due to its high pathogenicity and widespread occurrence around the world [1]. Owing to the lack of effective vaccines and other control measures, three classes of anthelmintics (benzimidazoles, macrocyclic lactones and cholinergic agonists) are dominant in this parasite control. Among these, ivermectin (IVM) has been successful in the control of many parasites in both humans and animals for more than 30 years [2,3]. However, as the uncontrolled use in livestock, the success is gradually undermined by the drug selection and spread of resistant parasite populations [4][5][6]. Consequently, its resistance has become widespread in many species of GIN of livestock.
Our knowledge about the resistant mechanisms of IVM has increased greatly over the last decades. Due to the technique limitations, the initial research on the IVM resistant mechanism mainly focused on the polymorphisms and expression level of candidate genes encoding the target receptors, such as Hc-lgc-37, Hc-glc-3, Hc-glc-5 and Hc-avr-14 [7][8][9]. Subsequently, attention turns to observing differences (polymorphisms and expression level) between resistant and susceptible strains or pre-and posttreatment in the P-gps and dyf-7 [10][11][12][13]. However, due to the extremely high levels of genetic diversity among H. contortus populations and lack of consistency between studies [14][15][16], the major genetic mediators of IVM resistance have not been unequivocally identi ed.
Fortunately, the high-quality reference genome of this parasite together with the increasing accessibility to high throughput RNA sequencing offers an opportunity to facilitate aimed studies to explore its resistant mechanism. Transcriptome sequencing has been used widely to identify candidate genes/pathways related to drug resistance [17,18]. However, the transcriptome changes of H. contortus IVM eld resistant strain induced by natural selection, have received little attention. In the present study, to explore transcriptomic difference and identify the candidate genes potentially related with drug resistance induced by long term selection of IVM in the eld sample, the H. contortus transcriptome of adult worms were sequenced and a comparative analysis was performed between IVM susceptible and eld resistant strains. The identi ed differentially expressed genes (DGEs) and subsequent enrichment analysis provided a theoretical basis for further understanding of IVM resistance in H. contortus. The control test and faecal egg count reduction test (FECRT) In order to further con rm the IVM resistance of Zhaosu strain, six goats (6-months old goats free of parasites) were infected with H. contortus of Zhaosu strain (7000 L3s per goat) and randomly divided into treated and untreated groups (three animals per group). The goats in treatment group were treated with 0.2 mg/kg (with sub-cutaneous injection) IVM thirty days after infection, and all animals were necropsied 14 days post-treatment and H. contortus worm burdens were counted. Resistance was con rmed when the reduction in mean worm counts was less than 90% (R (%) = (1-[E/C])×100 where E and C represent mean worm burden of treated and un-treated group) [19].
In addition, the FECRT was conducted based on the McMaster method. In brief, 2 g of individual faecal samples were homogenised in 58 ml of saturated sodium nitrate. Then, 0.15 ml of the suspension was added to each chamber of the slide for egg count after 5 min. The eggs per gram (EPG) was calculated by multiplying the number of eggs in the two chambers by 100. The FEC reduction was calculated according to the formula FECR (%) = 100 ×(1-[Xt/Xc]) where c and t mean the average of EPG for untreated and treated group, respectively, at 14 day post treatment, and resistance is con rmed when FECR is <95%, and the 95% con dence level is less than 90% [20].
Larval development assay (LDA) Fresh eggs were collected from the faeces based on a previous study [21]. LDA was carried out as modi ed from a previous study [22]. Brie y, the 2.4 ml suspension (2 ml egg suspension and 400 µl growth medium, ~5000 eggs) was added to T25 cell bottle and incubated at 27°C for 24 h for the hatching of the rst-stage larvae (L1). On the next day, the L1 larval suspension (99 µl, ~100 larvae) was aliquoted into the 96-well plates and 1 µl IVM working solution was added to each well and then the L1s were incubated at 27°C for another 6 days. IVM concentrations ranged from 12.5 ng/ml to 0.2 ng/ml for Zhaosu strain and 5 ng/ml to 0.1 ng/ml for Haecon-5 strain. The number of developed third-stage larvae (L3s) was counted and the L3 developmental rate was expressed as a percentage of the mean number in control wells. Three separate assays with triplicate were conducted for each concentration. The data were analysed using GraphPad Prism 8 software. Resistance ratio (RR) was calculated according to the formula: IC50 resistant isolate/IC50 susceptible isolate [23].

Sample preparation
The adult male and female worms of H. contortus were isolated from the abomasum of goats 30 days after infection with 7000 L3 of either Haecon-5 strain or Zhaosu strain. These worms were washed thoroughly in PBS and randomly assigned into four groups of 15 worms each with three biological replicates (resistance males, susceptible males, resistance females, susceptible females). All samples transferred to liquid nitrogen for storage until use.

RNA extraction, library preparation and sequencing
Total RNA was extracted using TRIzol® Reagent with manufacturer's instructions (Invitrogen, USA). RNA integrity was precisely assessed using 2100 Bio-analyzer (Agilent Technologies USA), and RNA samples with RIN≥8 were used to construct library. The transcriptome libraries were prepared using Illumina kit (San Diego, CA, USA) following the manufacturer's instructions. Brie y, polyA mRNA was isolated using oligo(dT) beads and cleaved in fragmentation buffer. Then, cDNA was synthesized using a SuperScript dscDNA synthesis kit (Invitrogen, CA) with random hexamer primers (Illumina). After quanti ed, the library was sequenced with the Illumina NovaSeq 6000 sequencer (Majorbio, Shanghai, China). The samples were named as RM1, RM2, RM3; RF1, RF2, RF3; SM1, SM2, SM3; SF1, SF2, SF3, respectively. The RNA-seq raw data were submitted to the National Centre for Biotechnology Information (NCBI) Sequence Read Archive with a BioProject ID: PRJNA772807.
To better evaluate the potential roles of the DEGs, the functional analysis of the DEGs was carried out based on the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG). GO and KEGG pathway enrichment analyses were performed via Majorbio cloud platform (https://cloud.majorbio.com/).

Validation of transcriptome sequencing results
Eleven DEGs were selected and primers designed using Primer-BLAST (Additional le 1: Table S1). Total RNA from the resistant adult male worms (RM), resistant adult female worms (RF), susceptible adult male worms (SM) and susceptible adult female worms (SF) was reverse-transcribed to the cDNA with HiScript II Q RT SuperMix (Vazyme, China). The qRT-PCR was conducted to assess mRNA levels of selected DEGs in biological triplicates with technical duplicates using a TB Green® Premix Ex Taq™ II (Takara, China). The reaction procedure: 95°C for 5 min, 40 cycles at 95°C for 15 s, 60°C for 30 s and 72°C for 15 s. The actin gene was used as an endogenous reference gene (Gene ID: DQ080917) [24], and relative expression fold change of each gene in resistant samples versus susceptible ones was calculated using the 2 −∆∆Ct method.

Statistical analysis
All experiments were independently conducted three times and analysed by GraphPad Prism 8.0.2 (GraphPad Company, USA). The t-test was performed to compare DEGs in resistant and susceptible strains and the results were de ned signi cant as P<0.05 (*P<0.05, **P<0.01, ***P < 0.001, ****P < 0.0001).

Detection of IVM resistance in H. contortus isolated from Zhaosu
The control test and faecal egg count reduction test (FECRT) The control and FECRT tests were carried out using 0.2 mg/kg of IVM in the treatment to assess IVM e cacy, respectively. The result of control test revealed that the percentage reduction in worm burdens was 32.5% (Table 1). In addition, the result of FECRT showed that the EPGs in treatment group were not decreased signi cantly at day 14 post treatment compared with the un-treatment group and the FECR was 1% in the IVM administration group (Table 2). C and E represent mean worm burden of un-treatment and treatment group. The results of LDA showed that with the increasing of the drug concentrations, larval developmental rate decreased and dropped to zero at the highest IVM concentration for each group (Fig. 1). The IC50 values were 0.218 ng/ml (95% CI: 0.208 to 0.227) and 1.291 ng/ml (95% CI: 1.209 to 1.377) for Haecon-5 and Zhaosu-R, respectively and the RR was calculated as 5.9.
Transcriptome sequencing and data assembly The total raw reads were generated from resistant and susceptible adult male and female worm transcriptomes (four groups) with three replicates for each sample, respectively. After quality control, a total of 93.98 Gb of clean data were obtained from the twelve samples altogether, each of which contained more than 7.19 Gb with quality score of Q20 ≥ 98.35 and Q30 ≥ 94.71, indicating high-quality sequencing (Additional le 1: Table S2). A range of 74.6-78.36% clean reads of each sample were aligned onto H. contortus reference genome (Additional le 1: Table S3). Principal component analysis (PCA) of the normalized RNA-seq read counts showed a high level of consistency and good separation between biological replicates of the same population from susceptible or resistant strains (Fig. 2) for hierarchical cluster analysis (Fig. 4a, b), showing that DGEs from the biological replicates were more closely clustered together. In addition, top ten up and down regulated genes in each comparison group were identi ed (Additional le 1: Table S4), containing many un-characterised genes (5 up and 7 down for RM, 2 up and 5 down for RF) and those with known or putative functions including parasitic stage speci c protein 1, RNA-directed DNA polymerase (Reverse transcriptase) domain containing protein, saposin type B domain containing protein, glutathione S-transferase, cytoplasmic dynein 2 light intermediate chain 1 and nematode cuticle collagen. In addition, among the DEGs identi ed in the transcriptomes, genes with functions involved in receptor activities, transport, detoxi cation, structural constituent of cuticle, lipid binding, transport, metabolism and signalling pathway were also identi ed.

DEGs encoding receptors
Early reports highlighted that IVM can interact with a wider range of ligand-gated channels found in parasitic nematodes, including GluClRs, γ-aminobutyric acid (GABA), nicotinic acetylcholine (nAChR) and glycine (Gly) receptors (Laing et al., 2017). Interestingly, among the DEGs identi ed in the present study, the GluClRs subunits encoding gene glc-5 (HCON_00161180) was down-regulated (log2 fold-change: -4.05) in RF. In addition, the genes that encoding the acetylcholine receptor were up and down regulated in the RM (up:2 down:4) and RF (up:1 down:1), and the G protein-coupled receptor (GPCR) genes were also up and down regulated in both RM (up:5 down:1) and RF (up:1 down:4) as compared with SM and SF (Additional le 1: Table S5).

DEGs involved in detoxi cation and transport
Except for the receptor encoding genes, 13 genes encoding molecules involved in detoxi cation including cytochromes P450 (CYPs), short-chain dehydrogenases/reductases (SDR), glutathione S-transferases (GSTs), and UDP-glycosyltransferases (UGTs), were up or down regulated (Additional le 1: Table S5). A total of four and two CYP450 encoding genes were found up-regulated in RM and RF, respectively. Two SDR genes were only up regulated in RM and three GST genes were down regulated both in RM and RF, respectively. One UGT gene was up regulated both in RM, another one and two UGT genes were up and down regulated in RF, respectively. In addition, the transcription level of the ABC transporter encoding gene P-gp-9.1 (HCON_00130050) was signi cantly higher in both RM and RF (log2 fold-change: 2.20 and 1.2, respectively), and Hco-abt-4 (HCON_00085890, log2 fold-change: 1.18) was signi cantly up-regulated in RF (Additional le 1: Table S5).

DEGs involved in lipid metabolism and cuticle collagen formation
Beyond mentioned above, most notable is that genes related with lipid metabolic process, cellular lipid metabolic process, lipid binding, fatty acid metabolic process, fatty acid biosynthetic process showed signi cantly up-regulated in RM (11 up-regulated) and RF (5 up-regulated) (Additional le: Table S6). In addition, 23 and 17 cuticle collagen genes showed differential expression in RM and RF, respectively, of which 21 (RM) and 15 (RF) were signi cantly down-regulated (Additional le 1: Table S6). In general, the GO term enrichment of the DEGs for the male and female showed a similar list of terms, especially in related with oxidation-reduction process, lipids and cuticle.

KEGG functional annotation of the DEGs
The ingenuity pathway analysis identi ed 79 (up-regulated DEGs) and 55 (down-regulated DEGs) for male, 96 (up-regulated DEGs) and 129 (down-regulated DEGs) pathways for female DEGs in resistant strain compared with susceptible strain SM and SF, respectively. The up-regulated DEGs related to KEGG pathway signi cantly enriched in autophagy-animal, lysosome, apoptosis, sphingolipid signalling pathway and ABC transporters in the male (Fig. 6a). The down-regulated DEGs signi cantly enriched in lysosome, autophagy-animal, apoptosis, antigen processing and presentation and ECM-receptor interaction (Fig. 6b).
In the female (Fig. 6), the up-regulated DEGs signi cantly enriched in nicotine addiction, central carbon metabolism in cancer, ABC transporters (Fig. 6c). The down-regulated DEGs signi cantly enriched in lysosome, apoptosis, sphingolipid signalling pathway, drug metabolism -cytochrome P450 (Fig. 6d). In general, the major enriched pathways for the male and female also showed a similar pattern in some pathway.  Fig. 7a, b).

Discussion
As H. contortus has shown a strong ability to develop resistance, IVM resistance in this parasite from different continents has been extensively studied and well documented [25][26][27][28][29][30][31]. However, our understanding of the molecular basis of the resistant mechanism is patchy. With advances in the application of genomic and transcriptomic technologies in exploring the mechanism of anthelmintic resistance, more recent studies focus on analysing the polymorphisms and expression level of the candidate genes. Therefore, in order to further understanding mechanism of IVM resistance in H. contortus, we carried out comparative transcriptome analysis between H. contortus adult worms of susceptible (Haecon-5) and resistant (Zhaosu) strains using RNA-Seq. In the present study, a total of 543 and 359 DEGs were identi ed in the RM vs SM (up: 338, down: 205) and RF vs SF (up: 105, down: 254), respectively. Notably, although both sexes have the same selection pressure, the amount of DEGs in the RM was signi cantly higher than in the RF, and more DEGs are up-regulated in the RM, but it is opposite in the SF. Considering IVM treatment causing the dysregulation of the genes expression related with embryo development in nematodes [32], these differences may be related with the effects of IVM on the eggs in utero of the female worms. Furthermore, these differences may also be related to the different mechanism of IVM resistance between male and female worms, and the reasons hopefully can be identi ed in the future when more data from eld resistant strains are available.
In previous studies, GluClRs gene families often considered to be involved in IVM resistance in H. contortus. For example, it has reported that haplotype frequency and transcription level changes of Hcoglc-5 may be related with IVM resistance [33,34]. Furthermore, the Hco-glc-5 was also down-regulated in female worms of both the MHco4(WRS) and MHco10(CAVR) resistant strains comparing with the IVMsusceptible strain [35]. Hence, it was expected that long term selection would lead to transcriptional changes of GluClR genes in the resistant strain. As expected, the transcription level of Hco-glc-5 was signi cantly reduced in RF. However, recent studies also found that there are no signi cant differences in loci of Hc-glc-5 between IVM-resistant and sensitive isolates and no evidence of introgression in either backcross [15,16]. Taken together, the inconsistent results from different studies may be caused by differences in the genetic background of different strains, and suggest that Hc-glc-5 may be related to IVM resistance in some strains, but more powerful evidence is required to de ne whether Hc-glc-5 to be major gene conferring IVM resistance.
Although IVM is often considered to only act on different GluCl channels of nematodes, intriguingly, there are growing evidences that it also acts on other ion channels [36]. For example, in the free-living nematode Caenorabditis elegans, recent studies demonstrated that IVM inhibits C. elegans muscle L-AChR receptors [37] and abamectin acts more complex antagonistic effect on nAChRs [38]. In H. contortus, a recent study identi ed 26 candidate genes in relation to H. contortus IVM resistance including nicotinic acetylcholine receptors (nAChRs) and G protein-coupled receptors (GPCRs) encoding genes [39]. In another study, several nAChRs and GPCRs genes were also up and down-regulated in both resistant strains [35]. In the present study, seven nAChRs genes and ten GPCR genes were differentially regulated in the RM and RF, respectively. Among these genes, Hco-acr-6 (HCON_0016650, predicted to be involved in cell communication, nervous system process, regulation of membrane potential and located in synapse.) was up-regulated in both RM and RF, which is in agreements with a previous report (Hco-acr-6 upregulated in IVM-resistant strain (MHco10)) [35]. This gene may be one of the important candidate genes given it is a member of ligand-gated ion channel gene families. Meanwhile, down-regulated GO terms include ligand-gated channel activity, G protein-coupled receptor signalling pathway, ion transport and ion transmembrane transporter activity, implying that IVM may also have a signi cant effect on nAChRs and GPCRs. However, due to the complexity of two receptor families, it is a considerable challenge for studying modes of IVM action on these receptors, and further functional validation is required.
In H. contortus, it has been reported that an up-regulation of ABC transporters had a direct correlation with IVM resistant phenotype [40] and the nal concentration of IVM in H. contortus [41]. In this study, P-gp-9.1 was signi cantly up-regulated in both RM and RF, and ABC transport pathway enriched in KEGG analysis (Fig. 6a, c). These results are consistent with a role for ABC transporters in IVM susceptibility in H. contortus, implying that ABC transporters (particularly P-gp) clearly contribute to IVM e ux in eld resistant strain. Together, these results raise the possibility that P-gp could be a potential component of a set of diagnostic tools for the detection of IVM-resistance in H. contortus. However, due to the multigenic nature of IVM resistance, more eld resistant strains are required to be sequenced to con rm the suitability.
H. contortus is armed with many metabolizing enzymes (42 CYPs, 44 SDR, 34 UGT and 28 GST) to protect against environmental toxins and anthelmintics [42]. The previous study showed that silencing the gene CYP (HCON_00143950) increased the sensitivity to IVM in larvae of H. contortus [39]. In the adult worms of this parasite, although IVM induced changes in the expression of CYPs [43], it can't be metabolized by xenobiotic-metabolizing enzymes (XME) [44]. In the present study, CYPs and SDR were up-regulated and GST was down-regulated in the resistant strain, whereas UGT was both up-and downregulated. At rst sight, these upregulated genes could be elucidated as being part of candidate genes that metabolize the IVM. Yet the UGTs and GSTs were downregulated which expected to be up-regulated in resistant strain. In C. elegans, RNA interference (RNAi) caused phenotypes of nine CYP genes associated with the embryonic lethality, uncoordinated movement, slow and stunted growth, and reduced fat storage [45]. Furthermore, CYPs were considered involved in fatty acid metabolism after exposure to IVM in the C. elegans IVM resistant strain [46]. Interestingly, in the present study, up-regulated XME genes were mainly involved in oxidation-reduction process, iron ion binding, heme binding, and up-regulated DEGs were mainly enriched in oxidation-reduction, lipid metabolic process, cellular lipid metabolic process, fatty acid metabolic process. These results suggest that the upregulated XME genes may play a role in detoxify endogenously accumulated metabolic compounds, such as lipid metabolism.
Lipids play indispensable roles in many aspects of intra-or inter-cellular signalling, cellular membranes, energy storage in the organisms [47]. The recent research reported that lipids modulate the activity and expression of e ux pumps, contributing to the development of the resistance in cancer [48,49]. In addition, the dramatic changes after exposure to IVM were dominated by genes related with lipid metabolism (C. elegans, resistant strain) and rapidly consumption of lipid stores (Globodera pallida) [46,50]. Intriguingly, in the present study, genes related with lipids and lipid metabolism underwent the greatest change in the up-regulated terms, include lipid metabolic process, fatty acid biosynthetic process, fatty acid metabolic process, lipase binding and activity, cellular lipid metabolic process, glycerolipid metabolic process terms (Fig. 5a, c). These great difference in energy homeostasis between the susceptible and resistant strains suggest that IVM has multiple effects involving both neuromuscular and non-neuromuscular targets. These up-regulated genes likely play key roles in nutritional requirements, anti-oxidation, membrane uidity, physical rigidity of the worm cuticle, to protect the worm itself from drug stimulation and adverse environment. To sum up, lipids and its metabolic pathways may play important roles in the IVM resistance, although it remains to be veri ed whether the lipids play a role in IVM-induced stress response, or direct role in IVM metabolism.
The cuticle is an external structure in all nematodes which plays a critical role in maintenance of body morphology and integrity, and locomotion [51]. Previous studies showed that the dysregulation of collagen encoding gene was considered to likely represent a non-speci c marker of stress, such as oxidation, response to multiple pathogens and drugs [51][52][53][54]. Furthermore, in C. elegans, it was found that trans-cuticular absorption may be a mode of entry for anti-nematodal IVM [55] and that dye-lling (dyf) mutations appear to affect cuticle permeability of IVM and the dye [56]. In addition, recent study reported that H. contortus can absorb nutrients and some anthelmintics (levamisole and macrocyclic lactones) through the cuticle [57]. In B. malayi, 17 out of 29 cuticle collagen genes were signi cantly down-regulated after exposure to IVM [32]. In the present study, 23 and 17 cuticle collagen genes displayed differential expression in both RM and RF, of which 21 and 14 were signi cantly downregulated (Table S6), respectively. In the GO enrichment analysis, collagen trimer and structural constituent of cuticle were signi cantly enriched in both RM and RF (Fig. 5b, d). These down-regulated genes play key roles the morphological changes in the integrity of the worm cuticle, and likely involved in protect the worm itself from IVM stimulation. These results suggest that the trans-cuticular penetration may be an important mode of IVM entry into worms and may represent novel candidate pathway for IVM resistant mechanism in the parasitic nematodes. The long-term IVM stimulation

Conclusions
In present study, we made comparative analyses on the transcriptome pro les between IVM resistant and susceptible strains of H. contortus. In addition to some candidate genes, we also found some DEGs related with nAChR, GPCRs, lipid metabolism and cuticle morphological formation. Our results suggested that there are different mechanisms in different populations and IVM resistance may be selected by changes of H. contortus metabolism. Our ndings provided new information for further research on the molecular mechanisms of the response to IVM in parasitic nematodes. Figure 1 Dose responses for the Haecon-5 and Zhaosu strain of Haemonchus contortus towards ivermectin. Data of each point from three separate experiments with triplicate (mean± SE).      Relative transcript levels of selected genes differentially expressed in male and female worms between IVM susceptible and resistant strains of Haemonchus contortus. a, b The validation of DEGs of male and worms between IVM susceptible and resistant strain. Data from the qRT-PCR represent the mean of three separate replicates, and bars represent standard error. *P < 0.05; **P < 0.01; *** P < 0.001, **** P < 0.0001.