Skip to main content
Download PDF

Ivermectin-induced gene expression changes in adult Parascaris univalens and Caenorhabditis elegans: a comparative approach to study anthelminthic metabolism and resistance in vitro

Parasites & Vectors volume 15, Article number: 158 (2022) Cite this article

Abstract

Background

The nematode Parascaris univalens is one of the most prevalent parasitic pathogens infecting horses but anthelmintic resistance undermines treatment approaches. The molecular mechanisms underlying drug activity and resistance remain poorly understood in this parasite since experimental in vitro models are lacking. The aim of this study was to evaluate the use of Caenorhabditis elegans as a model for P. univalens drug metabolism/resistance studies by a comparative gene expression approach after in vitro exposure to the anthelmintic drug ivermectin (IVM).

Methods

Twelve adult P. univalens worms in groups of three were exposed to ivermectin (IVM, 10–13 M, 10–11 M, 10–9 M) or left unexposed for 24 h at 37 °C, and total RNA, extracted from the anterior end of the worms, was sequenced using Illumina NovaSeq. Differentially expressed genes (DEGs) involved in metabolism, transportation, or gene expression with annotated Caernorhabditis elegans orthologues were identified as candidate genes to be involved in IVM metabolism/resistance. Similarly, groups of 300 adult C. elegans worms were exposed to IVM (10–9 M, 10–8 M and 10–7 M) or left unexposed for 4 h at 20 °C. Quantitative RT-PCR of RNA extracted from the C. elegans worm pools was used to compare against the expression of selected P. univalens candidate genes after drug treatment.

Results

After IVM exposure, 1085 DEGs were found in adult P. univalens worms but the relative gene expression changes were small and large variabilities were found between different worms. Fifteen of the DEGs were chosen for further characterization in C. elegans after comparative bioinformatics analyses. Candidate genes, including the putative drug target lgc-37, responded to IVM in P. univalens, but marginal to no responses were observed in C. elegans despite dose-dependent behavioral effects observed in C. elegans after IVM exposure. Thus, the overlap in IVM-induced gene expression in this small set of genes was minor in adult worms of the two nematode species.

Conclusion

This is the first time to our knowledge that a comparative gene expression approach has evaluated C. elegans as a model to understand IVM metabolism/resistance in P. univalens. Genes in P. univalens adults that responded to IVM treatment were identified. However, identifying conserved genes in P. univalens and C. elegans involved in IVM metabolism/resistance by comparing gene expression of candidate genes proved challenging. The approach appears promising but was limited by the number of genes studied (n = 15). Future studies comparing a larger number of genes between the two species may result in identification of additional candidate genes involved in drug metabolism and/or resistance.

Graphical Abstract

Introduction

The equine roundworm, Parascaris univalens, is one of the most important pathogenic parasites infecting foals, with a prevalence of 31–58% on stud farms [1,2,3,4]. Clinical manifestations of infection include stunted growth, respiratory symptoms, and in severe cases lethal obstruction and rupture of the small intestines [5]. In addition, P. univalens are known shedders of high numbers of eggs [6], and larvated eggs can remain viable on grazing fields for many years [7]. This asserts a high infection pressure on farms, resulting in foals requiring several treatments with deworming drugs during their first year [8].

Three main classes of anthelmintic drugs, benzimindzoles (BZ), tetrahydropyrimidines (THP), and macrocyclic lactones (MLs), are registered for use in horses to combat parasitic infections, including Parascaris spp. Macrocyclic lactones are the most commonly used drug class in veterinary medicine due to their high efficacy, low toxicity, and broad spectrum of target parasites [9]. The primary target of MLs in parasitic nematodes are glutamate-gated ion channels (GluCls), where for example the ML drug ivermectin (IVM) binds irreversibly to cause hyperpolarization of the pharynx and flaccid paralysis of muscles in the worm [10]. Development of anthelmintic resistance (AR) in Parascaris spp. has evolved as a result of extensive drug usage [11], and resistance to MLs is now widespread [3, 4, 11,12,13,14].

Several ML drug resistance mechanisms have been proposed, e.g. polymorphisms or changes in expression of the GluCls target and transport protein genes (see [15] for review). In addition, increased expression of genes encoding drug-metabolizing enzymes has been shown to be involved in resistance [16]. Polymorphisms cause conformational changes and alter drug-binding sites in the GluCl targets in the pathogenic nematodes Cooperia oncophora [17] and Haemonchus contortus [18]. In addition, decreased expression of drug target genes may lead to a reduction of drug-binding sites and thus reduced drug effectiveness [19, 20]. Transport protein genes, particularly encoding P-glycoproteins (Pgps) and other ATP-binding-cassette family members, have been implicated in anthelmintic drug resistance through the elimination of xenobiotic substances from the cell, thus preventing the drugs from reaching their target sites [21,22,23,24]. Multiple studies have also reported higher constitutive expressions of Pgps in resistant parasite strains compared to their susceptible counterparts [20, 21]. The metabolism of drugs is biphasic. In the first phase (Phase I), the drug is converted by oxidation, reduction, or hydrolysis to a more reactive compound that can be conjugated with an endogenous molecule such as glutathione or glucose in the second phase (Phase II). As a result, a soluble, inactive drug is generated that can be expelled from the cell [25]. Increased constitutive expression of Phase I [26] and Phase II metabolic genes [27] was reported in resistant H. contortus isolates after in vitro exposure with BZs. However, the exact roles of Phase I and II enzymes in development of ML resistance in parasitic nematodes remain unknown.

The genetic mechanisms underlying ML resistance in P. univalens are inadequately understood, and further research is hampered by the complex host-dependent life cycle. Moreover, in vivo studies using the parasite host would be more comprehensive, but are challenged by both ethical and financial constraints. Previous studies have only been conducted in an in vitro setting using live adult worms isolated from slaughtered horses or euthanized foals from research herds [28,29,30]. However, not only is this approach laborious and costly, but it is also inefficient because slaughter of foals is rare and the number of research herds is limited. In addition, the unnatural in vitro culture environment has been reported to have undesirable effect on gene expression of P. univalens [29]. Therefore, there is a need for the development of in vitro experimental models for P. univalens. Recently, the larval stage of P. univalens has been explored as a possible in vitro model, but authors reported variation in gene expression between larvae and adult P. univalens [31]. The free-living nematode Caenorhabditis elegans has been used as an in vitro model for parasitic nematodes, particularly those belonging to the same taxonomic clade, V, such as H. contortus [32,33,34,35]. Although P. univalens belongs to a different clade, III, previous studies have employed transgenic lines of C. elegans as models for interpretation of ML-resistance in P. univalens [28, 36, 37]. Factors that support C. elegans as a viable model for parasitic nematodes include possession of similar drug targets [38, 39] and being cheap and easy to maintain in the laboratory. In addition, C. elegans has a short life cycle of 3.5 days and a plethora of powerful genetic tools for manipulation are readily available [40]. To improve the likelihood of deriving conclusions from P. univalens research that uses C. elegans as a model organism, it is necessary to assess its applicability, at the very least in terms of gene expression following xenobiotic exposure.

In the current study, we investigated differentially expressed genes (DEGs), which encode drug targets, transporters, transcription regulators, or enzymes involved in drug metabolism in adult P. univalens after in vitro exposure to the ML substance, IVM. We further investigated the expression profile of the C. elegans orthologues of the above genes in adult C. elegans after in vitro exposure to IVM. Our objective was to evaluate whether C. elegans is a suitable model for P. univalens by studying and comparing gene expression of selected candidate genes in P. univalens and C. elegans after in vitro IVM exposure. Our findings showed that candidate genes, including the putative drug target gene, lgc-37, responded to IVM in P. univalens, but marginal to no response was observed in C. elegans despite dose-dependent behavioral effects observed in C. elegans after IVM exposure.

Materials and methods

Parascaris univalens

Adult P. univalens worms collected from two anthelmintic-naïve Icelandic foals at an abattoir in Selfoss, Iceland, were exposed to IVM in vitro. A detailed account of the experimental setup has been described in Martin et al. [29]. Care was taken to choose only the largest worms (i.e. females) for the study. Briefly, worms were incubated for 24 h in cell culture media containing 0.1% DMSO supplemented with 10–13, 10–11 or 10–9 M IVM. In addition, a control group (media+DMSO) of unexposed worms receiving cell culture media containing 0.1% DMSO was used (Fig. 1a). The experiment was performed in three biological replicates with at least three worms in each group. One worm from each subgroup was dissected. The anterior end (pharynx, nerve cell bodies (‘ganglia’), nerve ring, and part of the intestine) of the worm was cut into fine pieces and suspended in 1 ml Trizol (Invitrogen, Carlsbad, USA). After homogenization with a glass tissue grinder, chloroform was added and the aqueous phase of the homogenized suspension was isolated. The aqueous phase was advanced into NucleoSpin® RNA Plus Kit (Macherey Nagel, Düren, Germany) for RNA extraction according to the manufacturer’s instruction. Before preparing sequencing libraries from 500 ng total RNA, Fragment Analyzer (Agilent, Santa Clara, USA) was used to determine RNA concentration and quality. Illumina NovaSeq S1 flow cells and 100-bp paired end v1 sequencing chemicals were used to sequence three biological replicates per condition (Fig. 1a).

Fig. 1
figure 1

Illustration of the study design. Section a depicts experimental setup of Parascaris univalens. Adult worms were exposed to IVM drug concentrations, 10–13, 10–11, 10–9 M, and media+DMSO (control). After viability scoring, RNA was sequenced, candidate genes selected, and quantitative RT-qPCR performed on the selected candidate genes. Section b depicts experimental setup of Caenorhabditis elegans. Worm pools of adults at maximum reproduction were exposed to IVM drug concentrations, 10–9, 10–8, 10–7 M, and media+DMSO (control). Quantitative RT-PCR was performed on the P. univalens orthologues of selected candidate genes. Relative gene expression of the candidates was compared between nematode species

Full size image

RNA sequencing analysis

Reads were quality assessed using blastp [41] in a pipeline available at https://github.com/SLUBioinformaticsInfrastructure/RNAseq_nf. Mapping and quantification of reads against the predicted reference P. univalens transcriptome (https://parasite.wormbase.org/Parascaris_univalens_prjna386823/Info/Index/) were performed using Salmon v.0.11.3 [42]. Salmon output was matrix-summarized using R package tximport [43] for DESeq2 v.1.22.2 R package to determine DEGs as log2foldchange). Wald test derived P-values from DESeq2 were modified for multiple testing using the Benjamini-Hochberg approach-based p.adjust function in R base [44]. Differentially expressed genes were defined as those with an adjusted P-value < 0.05. Functional annotations of DEGs were identified by a BlastP search (e-value ≤ 10−5) of their respective protein sequences against the Swiss-Prot database.

To determine the number of genes shared by different drug concentrations, gene IDs were used to build comparative Venn diagrams in R v.3.6.1 using the venn.diagram function of the VennDiagram package [45]. The get.venn.partition function was used to retrieve gene IDs corresponding to each Venn diagram partition.

Principal component analysis

To analyze differences in the response to IVM exposure among individual worms, principal component analysis (PCA) was employed using scaled read counts > 0 for 10,830 genes from each worm (n = 12). The prcomp function in R was used on the scaled read count to perform PCA. To determine how much variation each principal component (PC) accounts for in the data, standard deviations generated by PCA were used to compute variances. The percentage variances were visualized in a scree plot using barplot function in R. The top two PCs, which account for most variation in the data, were used to build a PCA plot using ggplot2 R package.

Selection of candidate genes

The ten most DEGs in P. univalens were identified for each IVM concentration and considered as potential candidate genes. In addition, genes shared by two or more IVM concentrations were considered as potential candidates. Using the BiomaRt package in R, the potential candidate gene IDs were used to retrieve their respective ontologies and C. elegans orthologues from WormBase ParaSite (https://parasite.wormbase.org). From the pool of potential candidates, the selection of candidate genes was based on the presence of an orthologue gene in C. elegans and that they were functionally characterized as drug targets or involved in at least one of the following processes: xenobiotic metabolism, metabolite transportation, and gene expression regulation.

The similarity of P. univalens and corresponding C. elegans orthologues from WormBase were verified on the I-TASSER web-server [46] using their respective amino acid sequences.

Caenorhabditis elegans

Caenorhabditis elegans culture and maintenance

Wild-type N2 Bristol strain of C. elegans obtained from Caenorhabditis Genetics Center (CGC) was used in this study.

Worms were grown and maintained on OP50 Escherichia coli seeded Nematode Growth Media (NGM) at 20 °C as per standard methods [47].

Worm synchronization

Unless otherwise stated all M9 buffer [48] used was supplemented with 0.005% Tween20 to avoid worms/embryos adhesion to plastic ware.

Caenorhabditis elegans of different life stages were cultured on 10-cm-diameter NGM agar plates for 72 h to gravid adults. Worms were harvested with ice-cold M9 buffer into 15-ml Falcon tubes and synchronized by bleaching with 1.3% NaOCl and 0.5 M KOH, according to methods described by Porta-de-la-Riva et al. [49]. Embryo floatation was performed with 60% sucrose-0.1 M NaCl mixture (ratio: 1:1) and subsequently re-suspended in M9 buffer. Hatching into L1 larvae was done on NGM agar plates without bacteria for ~ 20 h to obtain a synchronous L1 population.

Approximately 860 synchronized L1 larvae per 10-cm-diameter OP50 E. coli seeded NGM plate were incubated at 20 °C for 76 h to an adult population at maximum reproduction stage (Fig. 1b).

Caenorhabditis elegans ivermectin exposure

Synchronized N2 adults at maximum reproduction (~ 300 worms/ml) were initially incubated in S-complete media supplemented with 10–9, 10–11, 10–13 M IVM (+ 0.025% DMSO final concentration), and media+DMSO (control). However, because the above-mentioned IVM concentrations had no discernible effect on the worms' behavior, the concentrations were increased to 10–7, 10–8, 10–9 M IVM (+ 0.025% DMSO final concentration). All treatments were set up in biological quadruplets supplemented with 1 mg/ml OP50 E. coli and incubated for 4 h at 20 °C (Fig. 1b).

Behavioral assays

After 4 h incubation at 20 °C, worms were phenotypically scored for each IVM concentration and control at room temperature. Scoring was based on three assays and performed with minor modifications as previously described by Johnson et al. [50].

The pharyngeal pumping assay was performed on three randomly selected worms per biological replicate. Worms were placed on 6-cm-diameter NGM plate containing a lawn of OP50 E. coli. Under an 80× magnification stereomicroscope (Nikon, Amstelveen, The Netherlands), the number of pharyngeal pumps per minute for each worm was counted.

The thrashing assay was performed on 12 worms, 3 randomly selected from a pool of at least 6 worms per biological replicate. Scoring was performed within their respective treatment solutions. Worms were video recorded for at least a minute. For scoring, a photographic image of the pool was taken at 0 s. Worms were randomly assigned numbers between one and n, where n is the total number of worms in the pool. Assigned numbers were randomized thrice, and three numbers were drawn. Worms with the drawn number were used for the scoring. The number of thrashings per minute was recorded. A thrash was defined as a complete sinusoidal motion from maximum to minimum amplitude and back.

In the dispersal assay, three to five worms from each treatment group were placed on the opposite edge of an OP50 lawn on a 6-cm-diameter NGM plate. Plates were scored after ~ 35 min based on the percentage of worms that left the inoculation spot and reached the OP50 lawn.

Caenorhabditis elegans RNA extraction, quality check, and cDNA synthesis

The remaining approximately 300 C. elegans worms were retrieved from the wells, washed twice in ice-cold M9 buffer, and centrifuged to pellet. The pellet was frozen in liquid nitrogen, ground with a pestle, and suspended in 900 μl Trizol (Invitrogen, Carlsbad, CA, USA) during continuous grinding. The aqueous phase of the homogenized suspension was isolated using chloroform and then advanced into the NucleoSpin® RNA Plus Kit (Macherey Nagel, Düren, Germany) for RNA extraction according to the manufacturer's instructions. RNA quality and quantity checks were performed on TapeStation 4150 (Agilent, Santa Clara, CA, USA) according to manufacturer’s instructions. Two micrograms of total RNA was reverse transcribed using SuperScript™ III First-Strand Synthesis System (Invitrogen, Carlsbad, CA, USA) to generate cDNA, which was then used in the qPCR assays.

Primer design and optimization

Coding sequences (CDS) for the candidate genes were retrieved from the P. univalens genome (accession no. PRJNA386823) in the WormBase ParaSite database [51]. Primers were designed on an open-source online application, primer3 version 0.4.0, available at http://bioinfo.ut.ee/primer3-0.4.0/, aiming for an amplicon size of approximately 500 bp, spanning at least two exons (Additional file 4: Table S1). For genes with multiple transcripts, only CDSs for transcripts with the highest read counts in the RNA sequencing (RNAseq) data were chosen. Coding sequences were aligned in CodonCode Aligner program (version 9.0.1, CodonCode Corporation) and unique or homologous regions selected for the design of primers. Polymerase chain reaction (PCR) was performed on a thermocycler (Applied Biosystems, Waltham, MA, USA) in a reaction mix consisting of 12.5 μl ToughMix (QuantaBio, Beverly, USA), 1 μl each 10-μM primer, 9.5 μl water, and 1 μl cDNA as a template, synthesized from 1 µg P. univalens RNA extract. The PCR program was as follows: initial denaturation at 95 °C for 5 min, 40 cycles of denaturation at 95 °C for 45 s, 55 °C annealing temperature for 45 s, elongation at 72 °C for 1 min, and final elongation at 72 °C for 5 min. Amplicons were submitted for Sanger sequencing at Macrogen (Amsterdam, The Netherlands). Resulting sequences were quality assessed and aligned to obtain consensus sequences in CodonCode Aligner. As a validation step, BLAST searches of consensus sequences against their corresponding CDS were performed using WormBase ParaSite. Thereafter, consensus sequences were used to design Reverse Transcriptase-qPCR (RT-qPCR) primers (Additional file 4: Table S1) using the same online platform as above, aiming for an amplicon size of 75–150 bp.

Similarly, CDS of C. elegans orthologues were retrieved from WormBase [52] using respective gene IDs. For orthologues with splice variants, primers were designed for CDS regions that were homologous to all variants. Likewise, amplicons were Sanger sequenced at Macrogen (Amsterdam, The Netherlands) and validated by BLAST searches against their corresponding CDS in the WormBase database. Reverse Transcriptase-qPCR assays were optimized for primer concentration, annealing temperature, and efficiency (Additional file 4: Table S1).

Relative gene expression of candidate genes using reverse transcriptase qPCR

Samples and non-template controls were run in duplicate reactions of 25 μl, composed of 12.5 μl Quantitect SYBR Green PCR mix (Qiagen, Hilden, Germany), 1 μl of each 10-μM primer, 8.5 μl water, and 2 μl (dilution factor: 5) cDNA as template on a CFX96 Touch PCR machine (Bio-Rad, Solna, Sweden). The RT-qPCR program was as follows: initial denaturation at 95 °C for 15 min, 40 cycles of denaturation at 95 °C for 15 s, assay specific annealing temperature for 30 s, and elongation at 72 °C for 30 s. For amplicon verification, a subsequent melt curve analysis was performed by temperature elevation from 60 to 95 °C, with increments of 0.5 °C for 5 s.

The Livak method [53] was used to determine the relative expression of the candidate genes. The gene expression was normalized to the geometric mean of two references genes for each species, i.e. actin, act-5 (PgR070_g023) and flyceraldehyde-3-phosphate dehydrogenase, gpd-1 (PgB20_g009) for P. univalens, and tba-1 (tubulin, alpha) and eif-3.C (eukaryotic initiation factor) for C. elegans and related to media+DMSO controls (Additional file 5: Table S2 and Additional file 6: Table S3). Gene expression was presented as log2 fold change. As a validation step, the log2 fold change values of candidate genes from RNAseq DEG analysis and RT-qPCR assays were imported into R and graphically compared in gglpot2’s scatterplot.

Statistical analysis

The average counts of behavioral assays were calculated in Microsoft Excel (Additional file 7: Table S4) and then imported into GraphPad Prism [Version 9.1.0 (221)] for statistical analysis and visualization. In a one-way ANOVA with Dunnett’s multiple comparisons test, average counts for the drug-treated worms were compared to media+DMSO controls and graphically represented.

Log2 fold changes of gene expression in the IVM exposed biological replicates were imported into GraphPad Prism [version 9.1.0 (221)] for each candidate gene and graphically represented. Graphical data are presented as mean values with a standard error of the mean (SEM) where applicable.

Results

Parascaris univalens RNAseq and differential gene expression

Adult P. univalens worms, collected from Icelandic horses, were treated with sublethal doses of IVM in vitro (10–13, 10–11, 10–9 M IVM for 24 h), and non-treated worms served as controls (see “Materials and methods”). Total RNA with RIN values ranging from 6.6 to 8.7 was sequenced on an Illumina Nova Seq. After read quality control with fastp, the number of reads per sample ranged between 63 and 132 million. Read sequences are available at https://www.ebi.ac.uk/ena [European Nucleotide Archive (ENA), accession no. PRJEB37010]. Read mapping against the reference transcriptome (https://parasite.wormbase.org/Parascaris_univalens_prjna386823/Info/Index/) yielded an average rate of 80–90% mapped reads, indicative of a high similarity between this transcriptomic data and reference transcriptome.

After PCA, 12 PCs were revealed of which PC1 and PC2 accounted for the most variation in the data, 35.5% and 14.1%, respectively (Fig. 2 and Additional file 1: Fig. S1). PC1 grouped worms into IVM treated and IVM untreated, except for worm 12. However, PC2 revealed variation in response within biological replicates but worms 1, 4 and 12 appeared away from their treatment groups (Fig. 2). Although effort was taken to select only females for the experiment, individual 2 in 10–11 M lacked an egg-containing uterus upon dissection. However, because this individual did not appear to be an outlier in the PCA plot (Fig. 2), it was included in the analysis.

Fig. 2
figure 2

Principal component analysis plot of Parascaris univalens based on scaled read counts of 10,830 genes showing variance among individual worms after exposure to three IVM drug concentrations (10−13 M, 10−11 M and 10−9 M) and media+DMSO (control). PC1 and PC2 represent the largest variances in the data, further portrayed by red and black dotted ellipses, respectively

Full size image

Following DESeq2 analysis of the 1085 genes (> 0.5-fold, adj. P-value < 0.05) differentially expressed across the three drug concentrations, 65% (703 genes) had known functions according to Swiss-Prot annotations. Twenty-one percent of the DEGs were found in 10–13 M, 17% in 10–11 M, and 62% in 10–9 M (Additional file 2: Fig. S2). The overall amplitude of gene expression changes after IVM treatment was small; 84 genes were upregulated more than twofold and 147 genes were downregulated more than twofold. There was also a quite high level of variability in gene expression between the different worms.

Selection of candidate genes for further characterization was based on the presence of an orthologue in C. elegans, level of differential gene expression and functional characterization as drug target or involvement in at least one of the following processes: xenobiotic metabolism, metabolite transportation, or gene expression regulation. Fifteen candidate genes were selected (11 up- and 4 downregulated) and further categorized into six groups: drug target, phase I metabolic enzyme, phase II metabolic enzyme, transporter, transcription regulator, and others (Table 1). The I-TASSER web server was used to verify P. univalens–C. elegans similarity for 80% and 87% of the candidate genes based on molecular function and biological process, respectively (Table 1).

Table 1 I-TASSER web server comparison of selected candidate genes’ orthologues in Caenorhabditis elegans and Parascaris univalens (geneIDs begin with “Pg”)
Full size table

Caenorhabditis elegans behavior assays after IVM exposure

Synchronized C. elegans N2 adults were initially exposed to similar IVM concentrations as adult P. univalens (10–9, 10–11, 10–13 M). However, as there were no observable phenotypic or gene expression effects on C. elegans with such low IVM concentrations (data not shown), the drug concentrations were increased to 10–9, 10–8, 10–7 M, respectively.

Increased IVM concentration showed an observable effect on the worms in all behavioral assays (Fig. 3). Ivermectin exposure also showed a significant dose-dependent reduction of the number of pharyngeal pumps per minute (Fig. 3a). In the thrashing assay, IVM concentrations (10–8 and 10–7 M) significantly reduced thrashes, with even immobile worms in the 10–7 M concentration demonstrated by the very low adjusted P-value (P ≤ 0.001) (Fig. 3b). The dispersal assay showed a similar dose-dependent trend, but was significant only in 10–7 M IVM (Fig. 3c). Worms incubated in media+DMSO (control) showed no observable effects and no significant behavioral differences (Fig. 3). The 10–7 M concentration was excluded from the qPCR assays, as we aimed for mechanisms at play only in sub-lethal conditions.

Fig. 3
figure 3

Behavioral assays for Caenorhabditis elegans. a Pharyngeal pumps per minute, b body thrashes per minute, and c dispersal after IVM exposure, 10–9 10–8, 10–7 M, and control (media+DMSO). Mean number (± standard error, SE) of viability scores in treatments were compared to those of the control in a one-way ANOVA (Dunnett’s multiple comparisons test) to elucidate significant differences

Full size image
Fig. 4
figure 4

Scatter plot comparing relative gene expression (log2 fold change) of Parascaris univalens candidate genes between two methods [ RNAseq (orange) and RT-qPCR (blue)] after exposure to 10–13, 10–11, and 10–9 M IVM for 24 h

Full size image

Expression of candidate genes in P. univalens and C. elegans after IVM exposure

Log2 fold change of candidate genes from RNAseq DEG analysis appeared consistent to those from RT-qPCR, with very few exceptions (Fig. 4). To verify differential gene expression in IVM-treated P. univalens and C. elegans, RNA with RIN values ranging from 6.6 to 9.8 were used in RT-qPCR assays, and all gene expressions described below are based on these assays. For simplicity, only C. elegans gene names are used for both P. univalens and C. elegans. Corresponding P. univalens gene names can be referenced in Table 11.

In general, although IVM appears to have an effect on the expression levels of candidate genes in P. univalens, this effect appears to vary among individual worms, agreeing with the observations seen in the RNAseq-based PCA plot (Fig. 2). In C. elegans, however, even if there is a strong phenotypic effect of IVM, it has marginal to no effect on the change in expression levels of any of the 15 selected candidate genes (Figs. 5, 6, 7, 8).

Fig. 5
figure 5

Relative gene expression (log2 fold change) of the GABA receptor subunit lgc-37 in a Parascaris univalens after exposure to 10–13, 10–11, and 10–9 M IVM for 24 h and b Caenorhabditis elegans after exposure to 10–9 and 10–8 M IVM for 4 h

Full size image
Fig. 6
figure 6

Relative gene expression (log2 fold change) of Phase I and II metabolic genes in Parascaris univalens (a, c) after exposure to 10–13, 10–11, and 10–9 M IVM for 24 h. Relative gene expression (log2 fold change) of Phase I and II metabolic genes in Caenorhabditis elegans (b, d) after exposure to 10–9 and 10–8 M IVM for 4 h

Full size image
Fig. 7
figure 7

Relative gene expression (log2 fold change) of transport protein genes in Parascaris univalens a after exposure to 10–13, 10–11, and 10–9 M IVM for 24 h. Relative gene expression (log2 fold change) of transport protein genes in Caenorhabditis elegans b after exposure to 10–9 and 10–8 M IVM for 4 h

Full size image
Fig. 8
figure 8

Relative gene expression (log2 fold change) of nhr-3 transcription regulator in Parascaris univalens a after exposure to 10–13, 10–11, and 10–9 M IVM for 24 h. Relative gene expression (log2 fold change) of nhr-3 transcription regulator in Caenorhabditis elegans b after exposure to 10–9 and 10–8 M IVM for 4 h

Full size image

The GABA receptor subunit gene, lgc-37, a putative drug target for IVM, showed a generally increased trend in change in expression in P. univalens after IVM exposure (Fig. 5a); however, increments in C. elegans were marginal (Fig. 5b).

The change in expression profiles of the Phase I metabolic genes in P. univalens showed non-uniform patterns. The short chain dehydrogenases (SDR) dhs-2 and dhs-27 showed an increased change in expression while F13D11.4 was decreased in P. univalens (Fig. 6a). The expression change of Phase II metabolic genes varied. The glucuronosyltransferase gene ugt-54 was decreased in P. univalens, while transferase C13A2.7 showed a trend of increased change in expression (Fig. 6c). However, in C. elegans, all genes generally appear not to be affected by the IVM exposure (Fig. 6b, d).

All investigated transport protein genes were enriched with major facilitator superfamily (MFS) members, except cup-4, which belongs to the ligand-gated ion channel family. Overall change in expression of transport protein genes showed a decreased trend in P. univalens and minor to no change in C. elegans (Fig. 7). Genes hmit-1.2/hmit-1.3 and Y17G128B.25 showed an increased change in expression in P. univalens (Fig. 7a).

Finally, in response to IVM, transcriptional regulator, nhr-3, showed an increased change in gene expression trend in P. univalens (Fig. 8a) but marginal to no change in expression was observed in C. elegans (Fig. 8b).

The neural signaling protein, pkc-1, showed an increased change in gene expression in P. univalens exposed to 10–11 and 10–9 M IVM (Additional file 3: Fig. S3A), but showed little to no change was observed in C. elegans (Additional file 3: Fig. S3B). Only the 10–13 M IVM showed marginal increased change in expression of ribulose-phosphate 3-epimerase gene F08F8.7 in P. univalens, and little to no change in C. elegans across all IVM concentration was observed in both species (Additional file 3: Fig. S3).

Discussion

Recent studies show that the equine roundworm P. univalens has developed resistance to the three major classes of anthelmintic drugs used and that the emergence of multidrug-resistant variants is apparent [4, 12, 54]. While the cause of drug resistance drug in parasites is poorly understood, the primary resistance mechanism to benzimidazoles has been better explored than those of other drug classes, particularly in clade V nematodes (see [15] for review). Consequently, most of the proposed resistance mechanisms have been derived from clade V nematodes and may therefore differ in Clade III nematodes such as P. univalens. It is critical to have a thorough understanding of P. univalens resistance mechanisms to potentiate species-specific interventions for resistance control and/or alleviation. This, however, is largely constrained by the complexity of its host-dependent life cycle and a lack of in vitro experimental models.

In this study, we evaluated the usage of C. elegans as a model for P. univalens by a comparative gene expression approach. Of the 1085 DEGs in P. univalens after IVM exposure, 231 DEGs had a fold change above twofold. This is however not unique to P. univalens as another clade III member, Brugia malayi, exposed to 10–7 M IVM for 24 h, demonstrated similar modest gene expression changes (14 genes upregulated and 8 genes downregulated more than twofold) [55]. While changes in gene expression greater than twofold are generally considered more important, expression changes of up to 0.5-fold have been proposed to be as equally meaningful in RNAseq assays if the false discovery rate is < 5% (adj. P-value 0.05), with at least three biological replicates used and data analyzed using the Deseq2 R program [56]. All DEGs in the current study followed the above recommendation and were therefore treated equally. Fifteen DEGs were selected as candidate genes and were further evaluated and compared with orthologous genes in C. elegans exposed to IVM.

Although IVM affected gene expression in P. univalens, there was little to no response in expression of the candidate genes in C. elegans, despite the significant dose-dependent behavioral effect of IVM on adult C. elegans worms. The premises for C. elegans as a model for research in parasitic nematodes have been thoroughly described in other studies [40, 57,58,59]. For example, Ascaris suum and C. elegans share 68.9% predicted genes [60], which supports C. elegans as a suitable model for the Ascarididae family. However, despite all the above, our results show differences in expression of the selected genes between C. elegans and Parascaris univalens after exposure to IVM. This observation can be attributed to a variety of factors. C. elegans and P. univalens are members of different clades, V and III, respectively, and they also have different life cycles and evolutionary histories [61]. As a result, genes and gene families associated with parasitism are either absent or have different functions in C. elegans [62,63,64]. The chosen candidate genes in this study are a small subset (1.4%) of all the genes responding to IVM exposure in P. univalens. Therefore, a global transcriptomic approach comparing whole worm transcriptomes of both C. elegans and P. univalens after IVM exposure may be a better method to evaluate whether C. elegans can be an accurate in vitro model for P. univalens.

The putative IVM target, lgc-37 (GABA receptor subunit), showed an increased expression in P. univalens, and Martin et al. [29] reported a similar observation. Increased expression of lgc-37 has not been correlated to AR; however, an amino acid substitution (K169R) in H. contortus lgc-37 transfected in C. elegans has been reported to exhibit reduced sensitivity to ML [65]. Taken together, these findings indicate that lgc-37 is important and that additional research is necessary to fully understand its role in P. univalens’ response to IVM.

In P. univalens, decreased expression was observed for the genes belonging to MFS-transport proteins slc-17.2 and F17C11.12 and the ligand-gated ion channel cup-4 after exposure to IVM. The role of MFS genes in parasitic worms has yet to be determined, but their involvement in drug resistance in bacteria and yeast has been reported [66]. Other than the efflux role, MFSs mediate cellular uptake of glucose and other saccharides [67]. Because IVM causes worm starvation through pharyngeal paralysis [68,69,70], we speculate that the decreased glucose concentration, a consequence of starvation, may trigger decreased gene expression of MFS responsible for glucose or other saccharide uptake. The ATP-binding cassette transporters, primarily Pgp efflux pumps, have been implicated in resistance in parasitic nematodes [21, 22, 24]. Interestingly, we found no DEGs of Pgps in our RNAseq dataset, which is consistent with previous studies in P. univalens [28, 30, 31, 71]. However, other studies have reported induced gene expression of Pgps in IVM-resistant larvae of Cooperia oncophora [72] and multi-drug resistant larvae of H. contortus exposed to IVM [22].

Differential expression of several Phase I enzymes was observed in P. univalens after IVM exposure. Expressions of dhs-2 and dhs-27 were increased whereas the reverse was seen for dhs-4 and F13D11.4. However, in C. elegans, all SDR genes showed a consistent minor decrease in expression levels. The role of SDRs in anthelmintic resistance is not well known. However, in several helminths, SDR genes have been shown to be involved in metabolic activities of BZs [73,74,75,76,77]. Even though we did not observe any major expressional response of SDR genes in C. elegans after IVM exposure, an increase in expression of dhs-23 in a BZ-resistant strain of C. elegans after in vitro exposure to BZ derivatives has been reported [35].

Expression of the Phase II enzyme, ugt-54, was decreased in P. univalens after IVM exposure, but no response was observed in C. elegans. Similarly, decreased expression of other members of the UGT family, ugt-3, ugt3a1, ugt-47, and ugt-48, have been observed in P. univalens in vitro exposed to pyrantel citrate, thiabendazole, and oxibendazole [29, 30]. Although we did not observe any response of ugt-54 in C. elegans, other studies have shown differential gene expression of UGTs in IVM-tolerant strain of C. elegans after exposure to 10–6 M IVM [78] and in the BZ-resistant strain of C. elegans after exposure to BZs [35, 79, 80]. Together, it is clear that UGTs are important in drug metabolism in both nematode species, and their role in anthelmintic resistance in P. univalens needs to be further elucidated.

The transcriptional regulator, nhr-3, showed increased change in expression in P. univalens after IVM exposure but a slightly decreased change in C. elegans. Increased expression of nhr-8 has reported to upregulate expression of Pgps and Phase I enzyme genes in C. elegans and resulted in a decreased efficacy of IVM [81]. A similar effect of IVM has also been observed in C. elegans carrying a transgenic construct expressing H. contortus-nhr-8 [81], suggesting that the function of nhr-8 could be similar in different nematode species. Even though nhr-3 does not respond in a similar way in C. elegans and P. univalens, the NHR family may play a role in AR.

Conclusion

In summary, this is the first time C. elegans has been evaluated as a model for P. univalens using a comparative gene expression approach. However, the genes selected showed dissimilar expression patterns between the two nematode species. This difference in expression could be due to caveats in the experimental design or simply the differences in life cycles of P. univalens and C. elegans. For example, due to the difficulties in obtaining live adult P. univalens, only three individuals were used per biological replicate in the study, compared to 300 worms for C. elegans. In addition, the unnatural in vitro culture environment may influence undesired behavioral changes or gene expression in P. univalens worms [29]. Furthermore, only 1.4% of RNAseq-derived DEGs from P. univalens were used as candidate genes in this study, of which some transcripts (splice variants) also showed varying expression patterns. Together, this limits the ability to draw firm conclusions regarding comparability between these two nematode species. However, with comparative genomics, using RNAseq data from IVM exposed C. elegans and P. univalens, it is plausible to find potential drug targets with similar expression in both nematodes. This could provide a better understanding of drug response mechanisms to further investigate and combat drug resistance.

Availability of data and materials

The pipeline used for quality assessment, mapping, and quantification is available (with nextflow and docker support) at https://github.com/SLUBioinformaticsInfrastructure/RNAseq_nf. Nucleotide sequence data reported in this paper are available in the European Nucleotide Archive (ENA) under the accession number PRJEB37010 (https://www.ebi.ac.uk/ena).

Abbreviations

IVM:

Ivermectin

BZ:

Benzimindzoles

THP:

Tetrahydropyrimidines

ML:

Macrocyclic lactones

GluCls:

Glutamate-gated ion channels

DEGs:

Differentially expressed genes

DMSO:

Dimethyl sulfoxide

PCA:

Principal component analysis

PC:

Principal component

ID:

Identification

NaOCl:

Sodium hypochlorite

KOH:

Potassium hydroxide

NaCl:

Sodium chloride

NGM:

Nematode growth media

CDS:

Coding sequences

PCR:

Polymerase chain reaction

RT-qPCR:

Reverse transcriptase quantitative polymerase chain reaction

ANOVA:

Analysis of variance

SEM:

Standard error of the mean

RIN:

RNA integrity number

GABA:

γ-Aminobutyric acid

SDR:

Short chain dehydrogenases

MFS:

Major facilitator superfamily

Pgp:

P-glycoprotein

References

  1. Laugier C, Sevin C, Ménard S, Maillard K. Prevalence of Parascaris equorum infection in foals on French stud farms and first report of ivermectin-resistant P. equorum populations in France. Vet Parasitol. 2012;188:185–9.

    Article  PubMed  Google Scholar 

  2. Relf VE, Morgan ER, Hodgkinson JE, Matthews JB. Helminth egg excretion with regard to age, gender and management practices on UK Thoroughbred studs. Parasitology. 2013;140:641–52.

    Article  CAS  PubMed  Google Scholar 

  3. Armstrong SK, Woodgate RG, Gough S, Heller J, Sangster NC, Hughes KJ. The efficacy of ivermectin, pyrantel and fenbendazole against Parascaris equorum infection in foals on farms in Australia. Vet Parasitol. 2014;205:575–80.

    Article  CAS  PubMed  Google Scholar 

  4. Alanazi AD, Mukbel RM, Alyousif MS, AlShehri ZS, Alanazi IO, Al-Mohammed HI. A field study on the anthelmintic resistance of Parascaris spp. in Arab foals in the Riyadh region, Saudi Arabia. Vet Q. 2017;37:200–5.

    Article  PubMed  Google Scholar 

  5. Cribb NC, Cote NM, Bouré LP, Peregrine AS. Acute small intestinal obstruction associated with Parascaris equorum infection in young horses: 25 cases (1985–2004). N Z Vet J. 2006;54:338–43.

    Article  CAS  PubMed  Google Scholar 

  6. Clayton HM, Duncan JL. The migration and development of Parascaris equorum in the horse. Int J Parasitol. 1979;9:285–92.

    Article  CAS  PubMed  Google Scholar 

  7. Clayton HM. Ascarids: recent advances. Vet Clin N Am Equine Pract. 1986;2:313–28.

    Article  CAS  Google Scholar 

  8. Nielsen MK. Evidence-based considerations for control of Parascaris spp. infections in horses. Equine Vet Educ. 2016;28:224–31.

    Article  Google Scholar 

  9. Shoop WL, Mrozik H, Fisher MH. Structure and activity of avermectins and milbemycins in animal health. Vet Parasitol. 1995;59:139–56.

    Article  CAS  PubMed  Google Scholar 

  10. Wolstenholme AJ, Rogers AT. Glutamate-gated chloride channels and the mode of action of the avermectin/milbemycin anthelmintics. Parasitology. 2005;131:S85–95.

    Article  CAS  PubMed  Google Scholar 

  11. Reinemeyer CR. Diagnosis and control of anthelmintic-resistant Parascaris equorum. Parasites Vectors. 2009;2:S8.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  12. Cooper LG, Caffe G, Cerutti J, Nielsen MK, Anziani OS. Reduced efficacy of ivermectin and moxidectin against Parascaris spp. in foals from Argentina. Vet Parasitol Reg Stud Rep. 2020;20: 100388.

    Google Scholar 

  13. Lindgren K, Ljungvall O, Nilsson O, Ljungström BL, Lindahl C, Höglund J. Parascaris equorum in foals and in their environment on a Swedish stud farm, with notes on treatment failure of ivermectin. Vet Parasitol. 2008;151:337–43.

    Article  CAS  PubMed  Google Scholar 

  14. Martin F, Svansson V, Eydal M, Oddsdóttir C, Ernback M, Persson I, et al. First report of resistance to ivermectin in Parascaris univalens in Iceland. J Parasitol. 2021;107:16–22.

    Article  PubMed  Google Scholar 

  15. Whittaker JH, Carlson SA, Jones DE, Brewer MT. Molecular mechanisms for anthelmintic resistance in strongyle nematode parasites of veterinary importance. J Vet Pharmacol Ther. 2017;40:105–15.

    Article  CAS  PubMed  Google Scholar 

  16. James CE, Hudson AL, Davey MW. Drug resistance mechanisms in helminths: is it survival of the fittest? Trends Parasitol. 2009;25:328–35.

    Article  CAS  PubMed  Google Scholar 

  17. Njue AI, Hayashi J, Kinne L, Feng XP, Prichard RK. Mutations in the extracellular domains of glutamate-gated chloride channel alpha3 and beta subunits from ivermectin-resistant Cooperia oncophora affect agonist sensitivity. J Neurochem. 2004;89:1137–47.

    Article  CAS  PubMed  Google Scholar 

  18. McCavera S, Rogers AT, Yates DM, Woods DJ, Wolstenholme AJ. An ivermectin-sensitive glutamate-gated chloride channel from the parasitic nematode Haemonchus contortus. Mol Pharmacol. 2009;75:1347–55.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. El-Abdellati A, De Graef J, Van Zeveren A, Donnan A, Skuce P, Walsh T, et al. Altered avr-14B gene transcription patterns in ivermectin-resistant isolates of the cattle parasites, Cooperia oncophora and Ostertagia ostertagi. Int J Parasitol. 2011;41:951–7.

    Article  CAS  PubMed  Google Scholar 

  20. Williamson SM, Storey B, Howell S, Harper KM, Kaplan RM, Wolstenholme AJ. Candidate anthelmintic resistance-associated gene expression and sequence polymorphisms in a triple-resistant field isolate of Haemonchus contortus. Mol Biochem Parasitol. 2011;180:99–105.

    Article  CAS  PubMed  Google Scholar 

  21. Dicker AJ, Nisbet AJ, Skuce PJ. Gene expression changes in a P-glycoprotein (Tci-pgp-9) putatively associated with ivermectin resistance in Teladorsagia circumcincta. Int J Parasitol. 2011;41:935–42.

    Article  CAS  PubMed  Google Scholar 

  22. Raza A, Kopp SR, Bagnall NH, Jabbar A, Kotze AC. Effects of in vitro exposure to ivermectin and levamisole on the expression patterns of ABC transporters in Haemonchus contortus larvae. Int J Parasitol Drugs Drug Resist. 2016;6:103–15.

    Article  PubMed  PubMed Central  Google Scholar 

  23. Schinkel AH, Jonker JW. Mammalian drug efflux transporters of the ATP binding cassette (ABC) family: an overview. Adv Drug Deliv Rev. 2003;55:3–29.

    Article  CAS  PubMed  Google Scholar 

  24. Xu M, Molento M, Blackhall W, Ribeiro P, Beech R, Prichard R. Ivermectin resistance in nematodes may be caused by alteration of P-glycoprotein homolog. Mol Biochem Parasitol. 1998;91:327–35.

    Article  CAS  PubMed  Google Scholar 

  25. Matoušková P, Vokřál I, Lamka J, Skálová L. The role of xenobiotic-metabolizing enzymes in anthelmintic deactivation and resistance in helminths. Trends Parasitol. 2016;32:481–91.

    Article  PubMed  CAS  Google Scholar 

  26. Yilmaz E, Ramünke S, Demeler J, Krücken J. Comparison of constitutive and thiabendazole-induced expression of five cytochrome P450 genes in fourth-stage larvae of Haemonchus contortus isolates with different drug susceptibility identifies one gene with high constitutive expression in a multi-resistant isolate. Int J Parasitol Drugs Drug Resist. 2017;7:362–9.

    Article  PubMed  PubMed Central  Google Scholar 

  27. Matoušková P, Lecová L, Laing R, Dimunová D, Vogel H, Raisová Stuchlíková L, et al. UDP-glycosyltransferase family in Haemonchus contortus: phylogenetic analysis, constitutive expression, sex-differences and resistance-related differences. Int J Parasitol Drugs Drug Resist. 2018;8:420–9.

    Article  PubMed  PubMed Central  Google Scholar 

  28. Janssen IJ, Krücken J, Demeler J, Basiaga M, Kornaś S, von Samson-Himmelstjerna G. Genetic variants and increased expression of Parascaris equorum P-glycoprotein-11 in populations with decreased ivermectin susceptibility. PLoS ONE. 2013;8: e61635.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Martin F, Dube F, Karlsson Lindsjö O, Eydal M, Höglund J, Bergström TF, et al. Transcriptional responses in Parascaris univalens after in vitro exposure to ivermectin, pyrantel citrate and thiabendazole. Parasites Vectors. 2020;13:342.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Scare JA, Dini P, Norris JK, Steuer AE, Scoggin K, Gravatte HS, et al. Ascarids exposed: a method for in vitro drug exposure and gene expression analysis of anthelmintic naïve Parascaris spp. Parasitology. 2020;147:659–66.

    Article  CAS  PubMed  Google Scholar 

  31. Martin F, Eydal M, Höglund J, Tydén E. Constitutive and differential expression of transport protein genes in Parascaris univalens larvae and adult tissues after in vitro exposure to anthelmintic drugs. Vet Parasitol. 2021;298: 109535.

    Article  CAS  PubMed  Google Scholar 

  32. Gibson SB, Harper CS, Lackner LL, Andersen EC. The Caenorhabditis elegans and Haemonchus contortus beta-tubulin genes cannot substitute for loss of the Saccharomyces cerevisiae beta-tubulin gene. MicroPubl Biol. 2021. https://doi.org/10.17912/micropub.biology.000411.

    Article  PubMed  PubMed Central  Google Scholar 

  33. Ondua M, Mfotie Njoya E, Abdalla MA, McGaw LJ. Investigation of anthelmintic activity of the acetone extract and constituents of Typha capensis against animal parasitic Haemonchus contortus and free-living Caenorhabditis elegans. Parasitol Res. 2021;120:3437–49.

    Article  PubMed  Google Scholar 

  34. Yates DM, Portillo V, Wolstenholme AJ. The avermectin receptors of Haemonchus contortus and Caenorhabditis elegans. Int J Parasitol. 2003;33:1183–93.

    Article  CAS  PubMed  Google Scholar 

  35. Stasiuk SJ, MacNevin G, Workentine ML, Gray D, Redman E, Bartley D, et al. Similarities and differences in the biotransformation and transcriptomic responses of Caenorhabditis elegans and Haemonchus contortus to five different benzimidazole drugs. Int J Parasitol Drugs Drug Resist. 2019;11:13–29.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Gerhard AP, Krücken J, Neveu C, Charvet CL, Harmache A, von Samson-Himmelstjerna G. Pharyngeal pumping and tissue-specific transgenic P-glycoprotein expression influence macrocyclic lactone susceptibility in Caenorhabditis elegans. Pharmaceuticals. 2021;14:153.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Janssen IJ, Krücken J, Demeler J, von Samson-Himmelstjerna G. Transgenically expressed Parascaris P-glycoprotein-11 can modulate ivermectin susceptibility in Caenorhabditis elegans. Int J Parasitol Drugs Drug Resist. 2015;5:44–7.

    Article  PubMed  PubMed Central  Google Scholar 

  38. Cully DF, Vassilatis DK, Liu KK, Paress PS, Van der Ploeg LH, Schaeffer JM, et al. Cloning of an avermectin-sensitive glutamate-gated chloride channel from Caenorhabditis elegans. Nature. 1994;371:707–11.

    Article  CAS  PubMed  Google Scholar 

  39. Dent JA, Smith MM, Vassilatis DK, Avery L. The genetics of ivermectin resistance in Caenorhabditis elegans. Proc Natl Acad Sci USA. 2000;97:2674–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Gilleard JS. The use of Caenorhabditis elegans in parasitic nematode research. Parasitology. 2004;128:S49–70.

    Article  CAS  PubMed  Google Scholar 

  41. Chen S, Zhou Y, Chen Y, Gu J. fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. 2018;34:i884–90.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  42. Patro R, Duggal G, Love MI, Irizarry RA, Kingsford C. Salmon provides fast and bias-aware quantification of transcript expression. Nat Methods. 2017;14:417–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Soneson C, Love MI, Robinson MD. Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences. F1000Res. 2015;4:1521.

    Article  PubMed  CAS  Google Scholar 

  44. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Ser B (Methodol). 1995;57:289–300.

    Google Scholar 

  45. Chen H, Boutros PC. VennDiagram: a package for the generation of highly-customizable Venn and Euler diagrams in R. BMC Bioinform. 2011;12:35.

    Article  Google Scholar 

  46. Yang J, Yan R, Roy A, Xu D, Poisson J, Zhang Y. The I-TASSER suite: protein structure and function prediction. Nat Methods. 2015;12:7–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Brenner S. The genetics of Caenorhabditis elegans. Genetics. 1974;77:71–94.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Stiernagle T. Maintenance of C. elegans. WormBook. 2006. p. 1–11.

  49. Porta-de-la-Riva M, Fontrodona L, Villanueva A, Cerón J. Basic Caenorhabditis elegans methods: synchronization and observation. J Vis Exp. 2012;64: e4019.

    Google Scholar 

  50. Johnson JR, Ferdek P, Lian LY, Barclay JW, Burgoyne RD, Morgan A. Binding of UNC-18 to the N-terminus of syntaxin is essential for neurotransmission in Caenorhabditis elegans. Biochem J. 2009;418:73–80.

    Article  CAS  PubMed  Google Scholar 

  51. Howe KL, Bolt BJ, Shafie M, Kersey P, Berriman M. WormBase ParaSite—a comprehensive resource for helminth genomics. Mol Biochem Parasitol. 2017;215:2–10.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Chen N, Harris TW, Antoshechkin I, Bastiani C, Bieri T, Blasiar D, et al. WormBase: a comprehensive data resource for Caenorhabditis biology and genomics. Nucleic Acids Res. 2005;33:D383–9.

    Article  CAS  PubMed  Google Scholar 

  53. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) method. Methods. 2001;25:402–8.

    Article  CAS  PubMed  Google Scholar 

  54. Martin F, Höglund J, Bergström TF, Karlsson Lindsjö O, Tydén E. Resistance to pyrantel embonate and efficacy of fenbendazole in Parascaris univalens on Swedish stud farms. Vet Parasitol. 2018;264:69–73.

    Article  CAS  PubMed  Google Scholar 

  55. Ballesteros C, Tritten L, O’Neill M, Burkman E, Zaky WI, Xia J, et al. The effects of ivermectin on Brugia malayi females in vitro: a transcriptomic approach. PLoS Negl Trop Dis. 2016;10: e0004929.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  56. Schurch NJ, Schofield P, Gierliński M, Cole C, Sherstnev A, Singh V, et al. How many biological replicates are needed in an RNA-seq experiment and which differential expression tool should you use? RNA. 2016;22:839–51.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  57. Hahnel SR, Dilks CM, Heisler I, Andersen EC, Kulke D. Caenorhabditis elegans in anthelmintic research—old model, new perspectives. Int J Parasitol Drugs Drug Resist. 2020;14:237–48.

    Article  PubMed  PubMed Central  Google Scholar 

  58. Holden-Dye L, Walker RJ. Anthelmintic drugs and nematicides: studies in Caenorhabditis elegans. WormBook. 2014. p. 1–29.

  59. Salinas G, Risi G. Caenorhabditis elegans: nature and nurture gift to nematode parasitologists. Parasitology. 2018;145:979–87.

    Article  PubMed  Google Scholar 

  60. Jex AR, Liu S, Li B, Young ND, Hall RS, Li Y, et al. Ascaris suum draft genome. Nature. 2011;479:529–33.

    Article  CAS  PubMed  Google Scholar 

  61. Blaxter ML, De Ley P, Garey JR, Liu LX, Scheldeman P, Vierstraete A, et al. A molecular evolutionary framework for the phylum Nematoda. Nature. 1998;392:71–5.

    Article  CAS  PubMed  Google Scholar 

  62. International Helminth Genomes Consortium. Comparative genomics of the major parasitic worms. Nat Genet. 2019;51:163–74.

    Article  CAS  Google Scholar 

  63. Viney M. The genomic basis of nematode parasitism. Brief Funct Genom. 2018;17:8–14.

    Article  CAS  Google Scholar 

  64. Wang J, Gao S, Mostovoy Y, Kang Y, Zagoskin M, Sun Y, et al. Comparative genome analysis of programmed DNA elimination in nematodes. Genome Res. 2017;27:2001–14.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  65. Feng XP, Hayashi J, Beech RN, Prichard RK. Study of the nematode putative GABA type-A receptor subunits: evidence for modulation by ivermectin. J Neurochem. 2002;83:870–8.

    Article  CAS  PubMed  Google Scholar 

  66. Paulsen IT, Brown MH, Skurray RA. Proton-dependent multidrug efflux systems. Microbiol Rev. 1996;60:575–608.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  67. Yan N. Structural advances for the major facilitator superfamily (MFS) transporters. Trends Biochem Sci. 2013;38:151–9.

    Article  CAS  PubMed  Google Scholar 

  68. Arena JP, Liu KK, Paress PS, Frazier EG, Cully DF, Mrozik H, et al. The mechanism of action of avermectins in Caenorhabditis elegans: correlation between activation of glutamate-sensitive chloride current, membrane binding, and biological activity. J Parasitol. 1995;81:286–94.

    Article  CAS  PubMed  Google Scholar 

  69. Avery L, Horvitz HR. Effects of starvation and neuroactive drugs on feeding in Caenorhabditis elegans. J Exp Zool. 1990;253:263–70.

    Article  CAS  PubMed  Google Scholar 

  70. Geary TG, Sims SM, Thomas EM, Vanover L, Davis JP, Winterrowd CA, et al. Haemonchus contortus: ivermectin-induced paralysis of the pharynx. Exp Parasitol. 1993;77:88–96.

    Article  CAS  PubMed  Google Scholar 

  71. Gerhard AP, Krücken J, Heitlinger E, Janssen IJI, Basiaga M, Kornaś S, et al. The P-glycoprotein repertoire of the equine parasitic nematode Parascaris univalens. Sci Rep. 2020;10:13586.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  72. De Graef J, Demeler J, Skuce P, Mitreva M, Von Samson-Himmelstjerna G, Vercruysse J, et al. Gene expression analysis of ABC transporters in a resistant Cooperia oncophora isolate following in vivo and in vitro exposure to macrocyclic lactones. Parasitology. 2013;140:499–508.

    Article  PubMed  CAS  Google Scholar 

  73. Cvilink V, Szotáková B, Krízová V, Lamka J, Skálová L. Phase I biotransformation of albendazole in lancet fluke (Dicrocoelium dendriticum). Res Vet Sci. 2009;86:49–55.

    Article  CAS  PubMed  Google Scholar 

  74. Cvilink V, Szotáková B, Vokrál I, Bártíková H, Lamka J, Skálová L. Liquid chromatography/mass spectrometric identification of benzimidazole anthelminthics metabolites formed ex vivo by Dicrocoelium dendriticum. Rapid Commun Mass Spectrom. 2009;23:2679–84.

    Article  CAS  PubMed  Google Scholar 

  75. Prchal L, Bártíková H, Bečanová A, Jirásko R, Vokřál I, Stuchlíková L, et al. Biotransformation of anthelmintics and the activity of drug-metabolizing enzymes in the tapeworm Moniezia expansa. Parasitology. 2015;142:648–59.

    Article  CAS  PubMed  Google Scholar 

  76. Solana HD, Rodriguez JA, Lanusse CE. Comparative metabolism of albendazole and albendazole sulphoxide by different helminth parasites. Parasitol Res. 2001;87:275–80.

    Article  CAS  PubMed  Google Scholar 

  77. Stuchlíková LR, Matoušková P, Vokřál I, Lamka J, Szotáková B, Sečkařová A, et al. Metabolism of albendazole, ricobendazole and flubendazole in Haemonchus contortus adults: sex differences, resistance-related differences and the identification of new metabolites. Int J Parasitol Drugs Drug Resist. 2018;8:50–8.

    Article  PubMed  PubMed Central  Google Scholar 

  78. Laing ST, Ivens A, Butler V, Ravikumar SP, Laing R, Woods DJ, et al. The transcriptional response of Caenorhabditis elegans to ivermectin exposure identifies novel genes involved in the response to reduced food intake. PLoS ONE. 2012;7: e31367.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  79. Jones LM, Rayson SJ, Flemming AJ, Urwin PE. Adaptive and specialised transcriptional responses to xenobiotic stress in Caenorhabditis elegans are regulated by nuclear hormone receptors. PLoS ONE. 2013;8: e69956.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  80. Laing ST, Ivens A, Laing R, Ravikumar S, Butler V, Woods DJ, et al. Characterization of the xenobiotic response of Caenorhabditis elegans to the anthelmintic drug albendazole and the identification of novel drug glucoside metabolites. Biochem J. 2010;432:505–14.

    Article  CAS  PubMed  Google Scholar 

  81. Ménez C, Alberich M, Courtot E, Guegnard F, Blanchard A, Aguilaniu H, et al. The transcription factor NHR-8: a new target to increase ivermectin efficacy in nematodes. PLoS Pathog. 2019;15: e1007598.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

Download references

Acknowledgements

The authors would like to express their gratitude to Sláturfélag Suurlands in Selfoss, Iceland, for allowing us to collect parasite material and to the Institute for Experimental Pathology at Keldur at the University of Iceland in Reykjavk, Iceland, for providing laboratory space. The SNP&SEQ Technology Platform in Uppsala performed the sequencing. The facility is a joint venture between Sweden’s National Genomics Infrastructure (NGI) and Science for Life Laboratory. Additionally, the Swedish Research Council and the Knut and Alice Wallenberg Foundation support the SNP&SEQ Platform. SLU Bioinformatics Infrastructure performed the bioinformatics analysis (SLUBI).

Funding

Open access funding provided by Swedish University of Agricultural Sciences. This work was supported by the Swedish Research Council FORMAS (Grant Number: 2018-01049).

Author information

Authors and Affiliations

  1. Division of Parasitology, Department of Biomedical Sciences and Veterinary Public Health, Swedish University of Agricultural Sciences, Box 7036, 750 07, Uppsala, Sweden

    Faruk Dube, Frida Martin & Eva Tydén

  2. Department of Cell and Molecular Biology, Uppsala University, 751 24, Uppsala, Sweden

    Andrea Hinas, Shweta Roy & Staffan Svärd

  3. Section of Immunology, Department of Biomedical Sciences and Veterinary Public Health, Swedish University of Agricultural Sciences, Box 7036, 750 07, Uppsala, Sweden

    Magnus Åbrink

Authors
  1. Faruk Dube
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Andrea Hinas
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Shweta Roy
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Frida Martin
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Magnus Åbrink
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Staffan Svärd
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Eva Tydén
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

FD, AH, MA, SS and ET designed the study. FD, FM, ET and SR collected parasite material. FD and SR performed the experiments. FD, ET AH, MA, and SS analyzed the data. FD wrote the manuscript. FM, AH, MA, SS, ET and SR reviewed the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Faruk Dube.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Figure S1.

A screeplot showing variation explained by all principal components.

Additional file 2: Figure S2.

Venn diagrams showing the number of differentially expressed genes in P. univalens that are shared among three IVM concentrations 10–13, 10–11, and 10–9 M. Genes with an adjusted P-value (Walds test and Benjamini-Hochberg procedure) < 0.05 were considered differentially expressed.

Additional file 3: Figure S3.

Relative gene expression (log2 fold change) of pkc-1 neural signaling protein and F08F8.7 Ribulose-phosphate-3-epimerase in Parascaris univalens (a) after exposure to 10–13, 10–11 and 10–9 M IVM for 24 h. Relative gene expression (log2 fold change) of pkc-1 neural signaling protein and F08F8.7 ribulose-phosphate-3-epimerase in Caenorhabditis elegans (b) after exposure to 10–9 and 10–8 M IVM for 4 h.

Additional file 4: Table S1.

Parascaris univalens and Caenorhabditis elegans primer sequences used in PCR assays as well as their respective amplicon sizes, annealing temperatures, and primer efficiencies.

Additional file 5: Table S2.

Parascaris univalens RT-qPCR cycle time values and log2 fold change calculations. IVM-13, IVM-11, IVM-9, and CD represent treatment conditions 10–13, 10–11, 10–9 M IVM, and media+DMSO (control), respectively.

Additional file 6: Table S3.

Caenorhabditis elegans RT-qPCR cycle time values and log2 fold change calculations. IVM-9, IVM-8, IVM-7, and IVM-D represent treatment conditions 10–9, 10–8, 10–7 M IVM, and media+DMSO (control), respectively.

Additional file 7: Table S4.

Behavioral assays of Caenorhabditis elegans after exposure to 10–9, 10–8, 10–7 M IVM and media+DMSO (control) for 4 h.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Dube, F., Hinas, A., Roy, S. et al. Ivermectin-induced gene expression changes in adult Parascaris univalens and Caenorhabditis elegans: a comparative approach to study anthelminthic metabolism and resistance in vitro. Parasites Vectors 15, 158 (2022). https://doi.org/10.1186/s13071-022-05260-4

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13071-022-05260-4

Keywords

Parasites & Vectors

ISSN: 1756-3305

Contact us