- Open Access
Transcriptional responses in Parascaris univalens after in vitro exposure to ivermectin, pyrantel citrate and thiabendazole
Parasites & Vectors volume 13, Article number: 342 (2020)
Parascaris univalens is a pathogenic parasite of foals and yearlings worldwide. In recent years, Parascaris spp. worms have developed resistance to several of the commonly used anthelmintics, though currently the mechanisms behind this development are unknown. The aim of this study was to investigate the transcriptional responses in adult P. univalens worms after in vitro exposure to different concentrations of three anthelmintic drugs, focusing on drug targets and drug metabolising pathways.
Adult worms were collected from the intestines of two foals at slaughter. The foals were naturally infected and had never been treated with anthelmintics. Worms were incubated in cell culture media containing different concentrations of either ivermectin (10−9 M, 10−11 M, 10−13 M), pyrantel citrate (10−6 M, 10−8 M, 10−10 M), thiabendazole (10−5 M, 10−7 M, 10−9 M) or without anthelmintics (control) at 37 °C for 24 h. After incubation, the viability of the worms was assessed and RNA extracted from the anterior region of 36 worms and sequenced on an Illumina NovaSeq 6000 system.
All worms were alive at the end of the incubation but showed varying degrees of viability depending on the drug and concentration used. Differential expression (Padj < 0.05 and log2 fold change ≥ 1 or ≤ − 1) analysis showed similarities and differences in the transcriptional response after exposure to the different drug classes. Candidate genes upregulated or downregulated in drug exposed worms include members of the phase I metabolic pathway short-chain dehydrogenase/reductase superfamily (SDR), flavin containing monooxygenase superfamily (FMO) and cytochrome P450-family (CYP), as well as members of the membrane transporters major facilitator superfamily (MFS) and solute carrier superfamily (SLC). Generally, different targets of the anthelmintics used were found to be upregulated and downregulated in an unspecific pattern after drug exposure, apart from the GABA receptor subunit lgc-37, which was upregulated only in worms exposed to 10−9 M of ivermectin.
To our knowledge, this is the first time the expression of lgc-37 and members of the FMO, SDR, MFS and SLC superfamilies have been described in P. univalens and future work should be focused on characterising these candidate genes to further explore their potential involvement in drug metabolism and anthelmintic resistance.
Nematodes within the genus Parascaris are pathogenic parasites of foals and yearlings worldwide. Traditionally the parasite has been referred to as Parascaris equorum but recent cytological studies have established that the major species currently infecting horses in the USA and Europe is Parascaris univalens [1,2,3]. Parascaris spp. infection causes nasal discharge, coughing and impaired growth, while large burdens can be lethal due to obstruction and perforation of the small intestine [4, 5]. To avoid parasite-related disease, most foals are usually treated with anthelmintics from the drug classes macrocyclic lactones, benzimidazoles or tetrahydropyrimidines several times during the first year . Macrocyclic lactones act by binding to parasite-specific glutamate- and γ-aminobutyric acid (GABA) gated ion channels in nerve and muscle cells, increasing the cells permeability to Cl- ions and leading to hyperpolarization which results in paralysis of the parasite . Benzimidazoles bind parasite β-tubulin molecules, thereby disrupting the polymerisation of microtubules, causing starvation and death of the worm . Tetrahydropyrimidines act as an agonist to the L-type nicotine acetylcholine gated ion channels, allowing Cl− to flow through, leading to depolarisation of muscle cells and spastic paralysis of the parasite .
Overuse of anthelmintic drugs has contributed to the development of resistance in several parasites of veterinary importance . The first reported case of anthelmintic resistance in Parascaris spp. was to the macrocyclic lactone ivermectin in 2002  and since then ivermectin resistance has been reported from around the world and is now considered widespread . Resistance to the tetrahydropyrimidine pyrantel was first discovered in the USA in 2008  and has also been found in Australia  and, more recently, in Europe [3, 15]. The benzimidazole fenbendazole is generally effective against Parascaris spp. in Europe [3, 15], but sporadic cases of treatment failure have been reported from Australia  and Saudi Arabia . Considering the risk of lethal complications in foals infected with Parascaris spp. and the lack of new anthelmintic drugs for the equine market, the development of resistance to all available drug classes in Parascaris spp. is a major threat to equine health and the equine industry.
Despite the increasing problem of anthelmintic resistance, the underlying causes of the development of resistance in parasitic nematodes are still poorly understood, particularly in ascarids. Advances in molecular techniques and genetics, such as the recent publication of a draft genome of P. univalens , have provided new possibilities to investigate responses in potential drug targets and pathways that may lead to resistance. There are several suggested mechanisms of anthelmintic resistance, including conformational changes or altered expression of the target molecule of the drug. Examples of conformational changes in the drug target are the three single nucleotide polymorphisms (SNPs) in the β-tubulin gene of strongyle nematodes. These mutations lead to amino acid substitutions and have been associated with the loss of action of benzimidazoles, particularly in strongyle nematodes of veterinary importance [18,19,20]. Decreased sensitivity to macrocyclic lactones in Haemonchus contortus was suggested to be connected to SNPs in the GABA receptor subunit gene lgc-37 [21, 22], although this was challenged . Additionally, reduced expression of the glutamate gated chloride channel avr-14 was described in ivermectin resistant isolates of Cooperia oncophora and Ostertagia ostertagi .
Other proposed mechanisms of anthelmintic resistance include changes in drug metabolising enzymes and alterations in efflux pumps leading to inactivation or removal of the drugs . Drug metabolism is usually divided into two phases; in phase I oxidation, reduction or hydrolysis convert the drug to a more reactive compound that can be conjugated with an endogenous molecule such as glutathione or glucose in phase II. This results in a soluble, inactive drug that can be removed from the cell . The phase I enzymes of the cytochrome P450-family (CYP) have been shown to be involved in drug resistance in insects , but the involvement of phase I and II enzymes in anthelmintic resistance of parasitic nematodes has not been thoroughly investigated, particularly not in ascarids. However, constitutively higher expression of a CYP34/35 family member  as well as the phase II enzyme uridine 5′-diphospho-glucuronosyltransferase-glucosyltransferase (UGT)  have been found in BZ resistant strains of H. contortus. Several studies have also shown evidence for upregulation of genes encoding drug efflux pumps, such as ATP-binding cassette (ABC) transporters, in macrocyclic lactone resistant strongyle nematodes [30,31,32]. Taken together, these results indicate that changes in drug metabolism and efflux, play a role in anthelmintic resistance.
Despite anthelmintic resistance in P. univalens being a growing threat to equine health, little is known about the molecular mechanisms behind the development of resistance. To our knowledge, no in vivo treatment experiments have been performed in this parasite to date and only two studies have reported the expression of genes after in vitro drug exposure of adult Parascaris spp. [33, 34].
The aim of our study was to identify genes that respond to drug treatment, focusing on drug targets, metabolising enzymes and transporters in adult P. univalens. We have used a whole transcriptome approach to compare the gene expression after in vitro exposure to sub-lethal doses of ivermectin, pyrantel citrate and thiabendazole.
Parasite material and karyotype
Adult P. univalens were collected from the intestines of two Icelandic foals, approximately six months-old at an abattoir in Selfoss, Iceland. The horses originated from the same farm in southern Iceland and had never been treated with anthelmintic drugs. After removal from the intestine, worms were rinsed with 37 °C PBS (Life Technologies, Carlsbad, USA) and transported to the laboratory (Institute for Experimental Pathology at Keldur, University of Iceland, Reykjavík, Iceland) in an insulated box. A faecal sample was taken from one of the foals for species identification by karyotyping of Parascaris spp. eggs as described in Martin et al. .
In vitro incubation experiment
To investigate the effect of in vitro incubation on viability and gene expression, two groups of worms were used. One group of nine worms was divided into three containers containing cell culture media (RPMI-1640 with the addition of 10% foetal bovine serum, 1% penicillin/streptomycin and 1% L-glutamine (Life Technologies, Carlsbad, USA) and then incubated for 24 h at 37 °C (control 24 h−DMSO). The other group of worms (n = 3) was immediately killed by removal of the anterior region upon arrival to the laboratory to serve as controls for the in vitro incubation (control 0 h) (Fig. 1).
In vitro drug exposure experiment
Worms were divided into three groups, one for each drug class. The worms were exposed to three different concentrations of ivermectin (IVM) (Sigma-Aldrich, Saint Louis, USA), pyrantel citrate (PYR) (Santa Cruz Biotechnology, Dallas, USA) or thiabendazole (TBZ) (Sigma-Aldrich), according to Table 1. The drug concentrations used were based on previous studies by Janssen et al.  for ivermectin exposure and Zhao et al.  for pyrantel citrate and thiabendazole exposure. All anthelmintic drugs were dissolved in DMSO (Swedish Veterinary Institute, Uppsala, Sweden) and serially diluted before addition to the media. For each concentration, three biological replicates, each consisting of three worms, were divided into different containers containing media and drug (final concentration of DMSO 0.1%) and incubated at 37 °C for 24 h. An additional group of three worms, serving as control, was incubated in media with the addition of 0.1% DMSO (control 24 h+DMSO) (Table 1, Fig. 1).
Viability scoring and dissection
After 24 h incubation, worm viability was assessed according to the scoring system developed by Scare et al.  to evaluate the effects of the different drug concentrations on viability. The worms in each container were observed and viability scores between 2 (movement only when stimulated with forceps) and 6 (seven or more whole body movements without stimulation) were given for each container. One worm from each container (n = 3) was dissected and the anterior region, containing the pharynx and a small part of the anterior intestine, was placed in individual tubes containing RNAlater (Ambion). Three worms were dissected directly upon arrival at the laboratory, serving as controls for the in vitro incubation experiment (control 0 h). The samples were then transported to Sweden (Swedish University of Agricultural Sciences, Department of Biomedical Sciences and Veterinary Public Health, Division of Parasitology) for the molecular analysis.
RNA extraction and sequencing
The samples were removed from RNAlater and cut into smaller pieces. One ml of Trizol (Invitrogen, Carlsbad, USA) was added to each sample and the tissue was homogenized in a glass tissue grinder. After addition of 0.2 ml of chloroform and centrifugation of the sample, 100 µl of the upper aqueous phase was mixed with 350 µl lysis buffer from the NucleoSpin® RNA Plus Kit (Macherey Nagel, Düren, Germany). RNA was then isolated according to the manufacturer’s instructions as described in the user manual of the kit. rDNase treatment and subsequent clean-up (NucleoSpin RNA Clean-up, Macherey Nagel) was performed according to the manufacturer’s protocol to ensure RNA purity.
Sample preparation and sequencing were performed at SciLifeLab Uppsala, SNP&SEQ Technology Platform. RNA concentration and integrity were checked by Fragment Analyzer (Agilent, Santa Clara, USA) before sequencing libraries were prepared from 500 ng total RNA using the TruSeq Stranded mRNA Library Preparation Kit including polyA selection (Illumina Inc., San Diego, USA) resulting in insert sizes of approximately 140 bp. Three biological replicates per condition were sequenced using Illumina NovaSeq S1 flow cells and 100 bp paired end v1 sequencing chemistry, resulting in 36 transcriptomes.
Read quality assessment, adaptor removal, filtering and removal of duplicates was performed using fastp  (Pipeline with nextflow and docker support available at https://github.com/SLUBioinformaticsInfrastructure/RNAseq_nf.). The resulting reads were mapped against the predicted transcriptome of P. univalens (parascaris_univalens.PRJNA386823.WBPS11.mRNA_transcripts.fa) available in WormBase ParaSite (https://parasite.wormbase.org/Parascaris_univalens_prjna386823/Info/Index/) and transcripts were quantified using Salmon v.0.11.3 . Transcript-level abundance, estimated counts and lengths from Salmon were summarised together with gene-ids from the P. univalens transcriptome into matrices by the R package tximport  for downstream differential expression analysis. Genes with five or more read counts were included in the differential expression analysis performed by the R package DESeq2, v.1.22.2 . P-values resulting from the Wald test incorporated into DESeq2 were adjusted for multiple testing using the Benjamini-Hochberg procedure  as applied by the R base p.adjust function. Genes with a log2 fold change of ≥ 1 or ≤ − 1 and an adjusted P-value < 0.05 were considered differentially expressed.
For principal components analysis (PCA), gene counts were transformed into log2 scale by applying regularized logarithm approach using the rlog function provided in the DESeq2 package and plotted using the plotPCA function. For all the above analyses, R v.3.5.2 was used . Functional annotations of differentially expressed genes were identified by searching the protein sequences against the Swiss-Prot database  using BlastP (e-value ≤ 10−5). Candidate genes for drug metabolism and drug targets significantly differentially expressed after exposure to one or more drugs were identified by their orthologue name. To identify the number of genes common among the concentrations of a particular drug, gene IDs were used to construct comparative Venn diagrams using the venn.diagram function of the VennDiagram package in R v.3.6.1 . Gene IDs corresponding to individual Venn diagram partitions were retrieved using the get.venn.partition function. Superfamilies of genes shared between concentrations were identified by BlastP searches (e-value ≤ 10−6) in NCBI (https://www.ncbi.nlm.nih.gov/).
Parascaris univalens was identified in the karyotype as DAPI stained eggs from the first mitotic division showed one pair of chromosomes.
Viability after in vitro incubation
After 24 h incubation, all worms were alive irrespective of drug and concentration used for exposure. Control worms were highly viable (score 6), whereas worms incubated in the highest drug concentrations were visibly less viable (score 3–4) than those incubated in lower concentrations (score 5–6) (Table 1).
Sequencing, quality control and mapping
RNA with RIN between 6.6–8.7 were sequenced on Illumina Nova Seq. The number of reads per sample varied between 63–132 million after quality assessment and filtering by fastp. Sequences have been deposited in the European Nucleotide Archive (ENA) under the accession number PRJEB37010 (https://www.ebi.ac.uk/ena). Mapping against the transcriptome resulted in an average mapping rate of 80–90% indicating high similarity between the transcriptomic data and the reference transcriptome.
Analysis of differentially expressed genes
The number of differentially expressed genes for each drug concentration and the proportion of uncharacterised genes are shown in Table 1. The PCA plot illustrates the differences in gene expression between individual worms (Fig. 2). It should be noted that the largest difference in gene expression was observed between the non-incubated worms (control 0 h) and both control and drug-exposed worms that were incubated in vitro for 24 h. Although care was taken to choose the largest worms (i.e. females) for the experiment, it was discovered during the dissection that two individuals (IVM 10−11 M individual 2 and TBZ 10−5 M individual 1) lacked egg-containing uterus. These two individuals did not show up as outliers in the PCA plot (Fig. 2) and were therefore included in the analysis. Log2 fold change and adjusted P-values for the differentially expressed genes discussed in this article are listed in Additional file 1: Table S1.
Transcriptional response in control worms after in vitro incubation
The largest number of differentially expressed genes was seen in the in vitro incubation experiment, where unexposed worms incubated in media (control 24 h−DMSO) were compared to non-incubated worms (control 0 h). Between these two groups, 1061 genes were differentially expressed and interestingly, among these were a number of candidate genes putatively involved in drug metabolism, such as CYPs, flavin containing monooxygenase (FMO), short chain dehydrogenase/reductase (SDR), glutathione S-transferase (GST) and UGTs, as well as ABC-transporters (Additional file 2: Table S2). Forty percent of the differentially expressed genes were uncharacterised.
Differential expression of putative drug targets
Putative drug targets were identified by their annotation and orthologue. Transcripts of ten drug targets were differentially expressed in worms exposed to IVM, PYR or TBZ compared to control 24 h+DMSO (Fig. 3a). An orthologue to the GABA-receptor subunit lgc-37 in H. contortus (GenBank: X73584), transcript PgR047_g061, was upregulated after exposure to the highest concentration of IVM (10−9 M). Three transcripts orthologous to acetylcholine receptor subunits (AChRs); PgB02X_g213, PgR034_g017 and PgR034_g018, were upregulated, whereas five AChRs; PgR043_g026, PgR075_g041, PgR005X_g207, PgR006_g054 and PgR045_g040, were downregulated according to Fig. 3a. The transcript PgE153_g002, identical to the previously described P. equorum β-tubulin isotype 2 (GenBank: KC713798), was downregulated after exposure to the highest dose of TBZ (10−5 M) but also after exposure to IVM (10−11 M) and PYR (10−8 M and 10−10 M) (Fig. 3a).
Comparison of the transcriptional response after exposure to different concentrations of IVM, PYR or TBZ
The ten most upregulated and downregulated genes after exposure to each drug are shown in Table 2. Of the most upregulated genes across all drugs, 57% are uncharacterized and only one transcript, the GABA-receptor subunit PgR047_ g061 (lgc-37), has previously been described as a drug target to IVM and to be involved in anthelmintic resistance. Of the most downregulated genes after exposure to the three different drugs, 53% are uncharacterized and two transcripts, the acetylcholine receptor PgR005X_g207 and the B-tubulin isotype 2 gene PgE153_g002, are previously described as drug targets of PYR and TBZ. No genes previously described to be involved in drug metabolism were found in this list.
Differentially expressed genes shared among the three concentrations of each drug are visualized in Fig. 4 and Table 3. After IVM exposure, 14 transcripts showed overlapping differential expression for all three concentrations (Fig. 4a). Of these, nine transcripts were upregulated whereas five were downregulated. The differentially expressed genes included four transcripts belonging to the SDR superfamily, putatively involved in phase I metabolism, two transcripts of the major facilitator superfamily (MFS), known to be involved in drug efflux and resistance in bacteria and yeast, and one transcript of the solute carrier superfamily (SCS), involved in transportation. Of the remaining seven differentially expressed transcripts, three do not have any known connection to drug metabolism or transport, while four transcripts are uncharacterised (Table 3).
After exposure to PYR, ten transcripts were commonly differentially expressed (Fig. 4b) and all but one were downregulated. As in the IVM exposed worms, the transcript PgR006_g137, a member of MFS, was downregulated (Table 3). Of the remaining differentially expressed transcripts, four belong to superfamilies not connected to drug metabolism or transportation, whereas five transcripts are uncharacterised (Table 3).
Six transcripts were commonly differentially expressed after exposure to all three concentrations of TBZ (Fig. 4c), five downregulated and one upregulated. Of these, two transcripts are involved in processes not connected to drug metabolism or transport, while four are uncharacterised or not similar to any known superfamily (Table 3).
Differential expression of other putative candidate genes
Putative candidate genes involved in drug metabolism were identified by their annotation and orthologue. Transcripts of 12 metabolising genes involved in phase I and phase II drug metabolism were identified in worms exposed to all drugs (Fig. 3b).
Four phase I transcripts were differentially expressed after drug exposure. One transcript belonging to the CYP family (PgR020_g037) was downregulated after exposure to the medium concentration of IVM (10−11 M). Three transcripts belonging to the flavin containing monooxygenase (FMO) superfamily were differentially expressed in an inconsistent pattern after drug exposure. The transcript PgR003_g196 was downregulated after exposure to IVM (10−9 M), TBZ (10−5 M) and PYR (10−10 M). The transcript PgR016_g066 was upregulated after exposure to IVM (10−9 M) only, whereas the transcript PgR161_g008 was upregulated after exposure to PYR (10−8 M) and TBZ (10−9 M) (Fig. 3b).
Eight genes belonging to two superfamilies of the phase II metabolizing pathway were also differentially expressed in an inconsistent pattern irrespective of drug and concentration used (Fig. 3b). One transcript, PgR024_g121, of the glutatione S-transferase superfamily (GST), was down regulated, but only after exposure to the lowest concentration of IVM (10−13 M). Seven transcripts of the uridine 5′-diphospho-glucuronosyltransferase superfamily (UGT) were downregulated in an unspecific pattern after exposure to all drugs (Fig. 3b).
The equine roundworm Parascaris spp. is resistant to several classes of anthelmintic drugs, but regardless of this, few studies focusing on the mechanisms involved have been published [33, 34, 45]. In this study, the transcriptomes of 36 individual P. univalens worms, confirmed by karyotyping, were analysed after in vitro exposure to different concentrations of IVM, PYR and TBZ. We found, to our knowledge for the first time in ascarid worms, that the expression of phase I gene families SDR and FMO were affected by exposure to anthelmintic drugs. We also found a 250-fold upregulation of a GABA receptor subunit.
The GABA receptor subunit PgR047_g061, orthologous to the H. contortus gene lgc-37, was upregulated after exposure to the highest IVM concentration (10−9 M) (Fig. 3). GABA receptors are targets for IVM in nematodes, but generally show lower affinity to the drug than the glutamate receptors . Even so, mutations in lgc-37 have been found to decrease the sensitivity of macrocyclic lactone drugs in H. contortus  suggesting a role in resistance. Transcriptional regulation of lgc-37 (PgR047_g061) in association with IVM response has not been reported previously, but taken together these results indicate that this gene could be an interesting candidate for further research in P. univalens.
Eight transcripts orthologous to AChR subunits, the putative drug target for PYR, were differentially expressed, but with no clear pattern after exposure to IVM, PYR and TBZ. Since the L-subtypes of AChRs have not yet been characterized in Parascaris spp., we cannot conclude whether these transcripts are the drug targets for pyrantel. Previous studies have shown that truncated forms of AChRs, as well as reduced expression of AChR subunits, have been associated with decreased sensitivity to PYR in strongyle nematodes [46, 47].
Furthermore, β-tubulin isotype 2 was downregulated after exposure to all three drug classes, however to a higher degree after TBZ exposure than after exposure to IVM or PYR. So far only two β-tubulin isotypes (isotype 1 and 2) have been described in Parascaris spp. , though more might be present in the genome since six putative β-tubulin genes were found in the closely related parasite Ascaridia galli . It has also been mentioned that the Ascaris suum genome contains at least nine β-tubulin genes . It should be noted that the isoforms of ascarid β-tubulins should not be directly compared to the strongyle nematodes as phylogenetic studies have shown that they are evolutionary separated [48, 49]. In a study by Martis et al. , where the transcriptomes of A. galli were compared before and after in vivo exposure to flubendazole, no differential expression of β-tubulins was found, whereas a study by Tyden et al.  showed an upregulation of β-tubulin isotype 1 in Parascaris spp. eggs exposed to thiabendazole during embryogenesis. However, downregulation of specific isotypes of β-tubulins has been observed in human cancer cells resistant to microtubule destabilising drugs [52, 53], indicating that expressional regulation of β-tubulin could be a cellular response to certain drugs.
Comparing the differential expression after exposure to three different drugs provided valuable insights on how expression of putative drug targets are similarly affected regardless of the drug used. Another point to consider is the inconsistent expression of putative drug targets and metabolising enzymes across the different concentrations of a drug. This could be explained by individual variation in gene expression between the biological replicates or can be an effect of changes in gene expression from drug metabolism to severe stress or apoptosis in response to high doses of anthelmintics. These are important findings to consider in the analysis of expression data and, to our knowledge, this is the first time this has been shown in P. univalens.
Several members of the phase I superfamily SDR were found to be differentially expressed after drug exposure in our study. The SDR family is a large superfamily present in all life forms , containing enzymes with conserved structure that metabolise endogenous and xenobiotic compounds by phase I reduction [55, 56]. Upregulation of SDR enzymes is believed to contribute to resistance towards chemotherapeutic agents in cancer treatment in human medicine  and has also been seen in Caenorhabditis elegans after in vitro exposure to benzimidazole drugs . The involvement of SDR in drug metabolism in parasitic nematodes is largely unknown. It has been shown that reduction is the main metabolic pathway for flubendazole metabolism in H. contortus, though it has not been established if the enzymes involved belonged to the SDR family or other related superfamilies . In summary, members of the SDR family may participate in the xenobiotic metabolism in helminths, but their possible involvement in anthelmintic resistance has not been studied in parasitic nematodes.
Members of the FMO family are phase I enzymes known to be involved in the biotransformation of xenobiotic compounds in many phyla . FMO enzymes have been shown to play a role in albendazole metabolism in the common liver fluke Fasciola hepatica and to be involved in the resistance to this drug [61, 62]. In contrast, Vokřál et al.  did not observe any activity of FMO enzymes in H. contortus in response to albendazole. In our study, two transcripts of the FMO family were upregulated after exposure to all drugs, indicating that they may play a universal role in the drug metabolism in P. univalens.
Only one member of the CYP family (PgR020_g037) was differentially expressed in this study, with a downregulation after exposure to IVM (10−11 M). This finding is in contrast to the upregulation of the CYP family member PgR071_g005 in P. univalens in vitro exposed to IVM and oxibendazole . Upregulation of CYP genes has also been observed in susceptible strains of H. contortus after in vitro exposure to IVM .
In summary, we found that several phase I enzymes were differentially expressed in P. univalens after drug exposure. To our knowledge, SDR and FMO enzymes have not yet been characterized in parasitic helminths and therefore their roles in drug metabolism or anthelmintic resistance have not been investigated. Although no consistent pattern was observed regardless of drug and concentration used, we found downregulation of several members of the phase II enzymes of the UGT family. In contrast, several studies have shown an increased expression and activity of UGTs in benzimidazole-resistant strains of H. contortus [29, 63]. Although the expression and function of phase I and phase II enzymes have been investigated in parasitic helminths, few studies have focused on ascarid worms and thus further investigation is required.
The upregulation of genes encoding transporters and drug efflux pumps have been suggested as a mechanism behind multidrug resistance in parasitic worms [30,31,32]. Focus has so far been mainly on certain ABC-transporters genes, but in agreement with Janssen et al. , we did not observe any changes in gene expression of pgp-11 and pgp-16 after drug exposure. On the other hand, we found that transcripts of the transportation superfamilies MFS and SCS were differentially expressed after exposure to anthelmintic drugs. These are large superfamilies of membrane proteins present in both prokaryotes and eukaryotes, transporting endogenous substances and drugs across cell membranes [64, 65]. The role of MFS proteins in drug resistance in parasitic helminths has so far not been investigated. However, their involvement in resistance to specific drugs, as well as multidrug resistance in bacteria and yeast, are well studied , indicating that these could be interesting candidates to investigate further in regards to anthelmintic resistance.
Among the ten most upregulated and downregulated genes after exposure to each of the three drugs, there were no previously described genes encoding xenobiotic enzymes or transporters. Even though the function of more than 50% of these genes is unknown, this list might contain novel candidates for drug metabolism and needs to be further investigated.
In a recently published study, Scare et al.  studied the transcriptomes of adult P. univalens after in vitro exposure to ivermectin and oxibendazole. Genes identified as potentially involved in drug detoxification and regulatory mechanisms differed from our results with only one similarity between the two studies, a member of the CYP family. These varying results are most likely due to the differences in study design between the two experiments, such as different drug concentrations and drug exposure times and the use of different tissues for RNA extraction. Scare et al.  also allowed a 24 h acclimatization period for the worms in tissue culture media before drug exposure, whereas we started the drug exposure at 0 h. These two studies show that despite similar aims, experimental design can affect differential expression considerably and thus also highlight the potential for future studies to identify important genes involved in drug interactions not found in either study.
In our experiment, unexposed control worms (control 24 h−DMSO) were not visibly affected by incubation at 37 °C for 24 h. Even so, we observed a large number (n = 1061) of differentially expressed genes in these worms when compared to non-incubated worms (control 0 h). Among these genes were several candidates putatively involved in drug metabolism and efflux. This implies that although not visibly affected by the in vitro incubation, the worms become stressed from the change in environment and nutritional supply when removed from the host intestine to incubation in media. This is an important point to consider when performing in vitro experiments. Although we cannot be sure that the in vitro response to drug treatment reflects what would happen in a similar in vivo situation, it is a way to obtain more knowledge about the parasite response to drug exposure. Since whole live animal experiments in foals are both expensive and pose ethical dilemmas, further development of in vitro models is crucial for research in parasitic nematodes.
We have explored the changes in gene expression across the complete transcriptome of adult P. univalens in response to in vitro exposure of the anthelmintic drugs ivermectin, pyrantel citrate and thiabendazole. We found a 250-fold upregulation of a possible target of ivermectin, an orthologue to the GABA receptor subunit lgc-37, which is an interesting candidate gene to investigate further. Surprisingly, few other candidate genes found in other parasitic worms, such as P-gps, were differentially expressed in our experiment. However, our data revealed differential expression of several novel gene candidates belonging to the phase I superfamilies SDR and FMO, as well as transporters belonging to the MFS and SLC families. To our knowledge, this is the first time the expression of these genes has been described in P. univalens and future work should be focused on characterising these candidate genes further to explore their potential involvement in drug metabolism and anthelmintic resistance.
Availability of data and materials
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).
short-chain dehydrogenase/reductase superfamily
major facilitator superfamily
solute carrier superfamily
single nucleotide polymorphisms
principal components analysis
RNA integrity number
acetylcholine receptor subunits
Nielsen MK, Wang J, Davis R, Bellaw JL, Lyons ET, Lear TL, et al. Parascaris univalens - a victim of large-scale misidentification? Parasitol Res. 2014;113:4485–90.
Jabbar A, Littlewood DTJ, Mohandas N, Briscoe AG, Foster PG, Muller F, et al. The mitochondrial genome of Parascaris univalens - implications for a “forgotten” parasite. Parasit Vectors. 2014;7:428.
Martin F, Hoglund J, Bergstrom TF, Karlsson Lindsjo O, Tyden E. Resistance to pyrantel embonate and efficacy of fenbendazole in Parascaris univalens on Swedish stud farms. Vet Parasitol. 2018;264:69–73.
Cribb NC, Cote NM, Boure 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.
Clayton HM, Duncan JL. Clinical signs associated with Parascaris equorum infection in worm-free pony foals and yearlings. Vet Parasitol. 1978;4:69–78.
National Veterinary Institute. Avmaskning mot spolmask. 2019. http://www.sva.se/djurhalsa/hast/parasiter-hos-hast/avmaskning-av-hast. Accessed 3 Mar 2020.
Martin R, Robertson A, Wolstenholme A. Mode of action of the macrocyclic lactones. In: Vercruysse J, Rew RS, editors. Macrocyclic lactones in antiparasitic therapy. Wallingford: CABI Publishing; 2002. p. 125–40.
Lacey E. Mode of action of benzimidazoles. Parasitol Today. 1990;6:112–5.
Martin RJ, Robertson AP. Mode of action of levamisole and pyrantel, anthelmintic resistance, E153 and Q57. Parasitology. 2007;134:1093–104.
Reinemeyer CR. Diagnosis and control of anthelmintic-resistant Parascaris equorum. Parasit Vectors. 2009;2(Suppl. 2):S8.
Boersema JH, Eysker M, Nas JWM. Apparent resistance of Parascaris equorum to macrocyclic lactones. Vet Rec. 2002;150:279–81.
Peregrine AS, Molento MB, Kaplan RM, Nielsen MK. Anthelmintic resistance in important parasites of horses: does it really matter? Vet Parasitol. 2014;201:1–8.
Lyons ET, Tolliver SC, Ionita M, Collins SS. Evaluation of parasitical activity of fenbendazole, ivermectin, oxibendazole, and pyrantel pamoate in horse foals with emphasis on ascarids (Parascaris equorum) in field studies on five farms in central Kentucky in 2007. Parasitol Res. 2008;103:287–91.
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.
Hautala K, Nareaho A, Kauppinen O, Nielsen MK, Sukura A, Rajala-Schultz PJ. Risk factors for equine intestinal parasite infections and reduced efficacy of pyrantel embonate against Parascaris sp. Vet Parasitol. 2019;273:52–9.
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.
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.
Kwa MS, Veenstra JG, Roos MH. Benzimidazole resistance in Haemonchus contortus is correlated with a conserved mutation at amino acid 200 in beta-tubulin isotype 1. Mol Biochem Parasitol. 1994;63:299–303.
Silvestre A, Cabaret J. Mutation in position 167 of isotype 1 beta-tubulin gene of trichostrongylid nematodes: role in benzimidazole resistance? Mol Biochem Parasitol. 2002;120:297–300.
Ghisi M, Kaminsky R, Maser P. Phenotyping and genotyping of Haemonchus contortus isolates reveals a new putative candidate mutation for benzimidazole resistance in nematodes. Vet Parasitol. 2007;144:313–20.
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.
Beech R, Levitt N, Cambos M, Zhou S, Forrester SG. Association of ion-channel genotype and macrocyclic lactone sensitivity traits in Haemonchus contortus. Mol Biochem Parasitol. 2010;171:74–80.
Laing R, Maitland K, Lecová L, Skuce PJ, Tait A, Devaney E. Analysis of putative resistance gene loci in UK field populations of Haemonchus contortus after 6 years of macrocyclic lactone use. Int J Parasitol. 2016;46:621–30.
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.
James CE, Hudson AL, Davey MW. Drug resistance mechanisms in helminths: is it survival of the fittest? Trends Parasitol. 2009;25:328–35.
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.
Daborn P, Boundy S, Yen J, Pittendrigh B, Ffrench-Constant R. DDT resistance in Drosophila correlates with Cyp6g1 over-expression and confers cross-resistance to the neonicotinoid imidacloprid. Mol Genet Genomics. 2001;266:556–63.
Yilmaz E, Ramunke S, Demeler J, Krucken 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.
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.
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.
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.
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.
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 naive Parascaris spp. Parasitology. 2020;147:659–66.
Janssen IJ, Krucken J, Demeler J, Basiaga M, Kornas 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.
Zhao J, Williams AR, Hansen TVA, Thamsborg SM, Cai J, Song S, et al. An in vitro larval migration assay for assessing anthelmintic activity of different drug classes against Ascaris suum. Vet Parasitol. 2017;238:43–8.
Scare JA, Steuer AE, Shaffer CL, Slusarewicz P, Mousley A, Nielsen MK. Long live the worms: methods for maintaining and assessing the viability of intestinal stages of Parascaris spp. in vitro. Parasitology. 2018;146:1–9.
Chen S, Zhou Y, Chen Y, Gu J. fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. 2018;34:i884–90.
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.
Soneson C, Love MI, Robinson MD. Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences. F1000Res. 2015;4:1521.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.
Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Series B Stat Methodol. 1995;57:289–300.
R Development Core Team. R: a language and environment for statistical computing. Vienna: R Foundation for Statistical Computing; 2018. https://www.R-project.org/.
UniProt Consortium. UniProt: a worldwide hub of protein knowledge. Nucleic Acids Res. 2019;47:D506–15.
Chen H, Boutros PC. VennDiagram: a package for the generation of highly-customizable Venn and Euler diagrams in R. BMC Bioinform. 2011;12:35.
Tyden E, Dahlberg J, Karlberg O, Hoglund J. Deep amplicon sequencing of preselected isolates of Parascaris equorum in beta-tubulin codons associated with benzimidazole resistance in other nematodes. Parasit Vectors. 2014;7:410.
Kopp SR, Coleman GT, Traub RJ, McCarthy JS, Kotze AC. Acetylcholine receptor subunit genes from Ancylostoma caninum: altered transcription patterns associated with pyrantel resistance. Int J Parasitol. 2009;39:435–41.
Neveu C, Charvet CL, Fauvin A, Cortet J, Beech RN, Cabaret J. Genetic diversity of levamisole receptor subunits in parasitic nematode species and abbreviated transcripts associated with resistance. Pharmacogenet Genomics. 2010;20:414–25.
Tyden E, Engstrom A, Morrison DA, Hoglund J. Sequencing of the beta-tubulin genes in the ascarid nematodes Parascaris equorum and Ascaridia galli. Mol Biochem Parasitol. 2013;190:38–43.
Martis MM, Tarbiat B, Tyden E, Jansson DS, Hoglund J. RNA-Seq de novo assembly and differential transcriptome analysis of the nematode Ascaridia galli in relation to in vivo exposure to flubendazole. PLoS One. 2017;12:e0185182.
Kotze AC, Hunt PW, Skuce P, von Samson-Himmelstjerna G, Martin RJ, Sager H, et al. Recent advances in candidate-gene and whole-genome approaches to the discovery of anthelmintic resistance markers and the description of drug/receptor interactions. Int J Parasitol Drugs Drug Resist. 2014;4:164–84.
Tyden E, Skarin M, Andersson-Franko M, Sjoblom M, Hoglund J. Differential expression of beta-tubulin isotypes in different life stages of Parascaris sppafter exposure to thiabendazole. Mol Biochem Parasitol. 2016;205:22–8.
Cheung CH, Wu SY, Lee TR, Chang CY, Wu JS, Hsieh HP, et al. Cancer cells acquire mitotic drug resistance properties through beta I-tubulin mutations and alterations in the expression of beta-tubulin isotypes. PLoS One. 2010;5:e12564.
Kavallaris M, Kuo DY, Burkhart CA, Regl DL, Norris MD, Haber M, et al. Taxol-resistant epithelial ovarian tumors are associated with altered expression of specific beta-tubulin isotypes. J Clin Invest. 1997;100:1282–93.
Kallberg Y, Oppermann U, Persson B. Classification of the short-chain dehydrogenase/reductase superfamily using hidden Markov models. FEBS J. 2010;277:2375–86.
Matsunaga T, Shintani S, Hara A. Multiplicity of mammalian reductases for xenobiotic carbonyl compounds. Drug Metab Pharmacokinet. 2006;21:1–18.
Kisiela M, El-Hawari Y, Martin HJ, Maser E. Bioinformatic and biochemical characterization of DCXR and DHRS2/4 from Caenorhabditis elegans. Chem Biol Interact. 2011;191:75–82.
Soldan M, Netter KJ, Maser E. Induction of daunorubicin carbonyl reducing enzymes by daunorubicin in sensitive and resistant pancreas carcinoma cells. Biochem Pharmacol. 1996;51:117–23.
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.
Cvilink V, Kubícek V, Nobilis M, Krízová V, Szotáková B, Lamka J, et al. Biotransformation of flubendazole and selected model xenobiotics in Haemonchus contortus. Vet Parasitol. 2008;151:242–8.
Cashman JR. Monoamine oxidases and flavin-containing monooxygenases. In: McQueen CA, editor. Comprehensice toxicology. Amsterdam: Elsevier; 2018. p. 87–125.
Alvarez LI, Solana HD, Mottier ML, Virkel GL, Fairweather I, Lanusse CE. Altered drug influx/efflux and enhanced metabolic activity in triclabendazole-resistant liver flukes. Parasitology. 2005;131:501–10.
Brennan GP, Fairweather I, Trudgett A, Hoey E, McCoy, McConville M, et al. Understanding triclabendazole resistance. Exp Mol Pathol. 2007;82:104–9.
Vokřál I, Jirásko R, Stuchlíková L, Bártíková H, Szotáková B, Lamka J, et al. Biotransformation of albendazole and activities of selected detoxification enzymes in Haemonchus contortus strains susceptible and resistant to anthelmintics. Vet Parasitol. 2013;196:373–81.
Paulsen IT, Brown MH, Skurray RA. Proton-dependent multidrug efflux systems. Microbiol Rev. 1996;60:575–608.
Hoglund PJ, Nordstrom KJ, Schioth HB, Fredriksson R. The solute carrier families have a remarkably long evolutionary history with the majority of the human families present before divergence of Bilaterian species. Mol Biol Evol. 2011;28:1531–41.
Open access funding provided by Swedish University of Agricultural Sciences. The authors are grateful to Sláturfélag Suðurlands, Selfoss, Iceland for allowing us to collect parasite material and Institute for Experimental Pathology at Keldur, University of Iceland, Reykjavík, Iceland for providing laboratory space. Sequencing was performed by the SNP&SEQ Technology Platform in Uppsala. The facility is part of the National Genomics Infrastructure (NGI) Sweden and Science for Life Laboratory. The SNP&SEQ Platform is also supported by the Swedish Research Council and the Knut and Alice Wallenberg Foundation. Bioinformatics analysis was performed by SLU Bioinformatics Infrastructure (SLUBI). We would also like to thank Dr Tom Martin for proof reading the manuscript.
This work was supported by the Swedish research council FORMAS (grant no. 942-2015-508).
Ethics approval and consent to participate
No ethical permissions were necessary for this study as the parasites were collected from horses slaughtered for meat production.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Additional file 1: Table S1.
Log2 fold change and adjusted P-values for the differentially expressed genes discussed in this article.
Additional file 2: Table S2.
Parascaris univalens in vitro incubation experiment. Differentially expressed genes (fold change ≥ 2 or and adjusted P-value (Walds test and Benjamini-Hochberg procedure) < 0.05) and their corresponding orthologues according to Swiss-Prot database after 24 h in vitro incubation in media (control 24 h-DMSO). Gene expression was normalized to non-incubated worms (control 0 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.
About this article
Cite this article
Martin, F., Dube, F., Karlsson Lindsjö, O. et al. Transcriptional responses in Parascaris univalens after in vitro exposure to ivermectin, pyrantel citrate and thiabendazole. Parasites Vectors 13, 342 (2020). https://doi.org/10.1186/s13071-020-04212-0
- Anthelmintic resistance
- RNA sequencing
- Differential expression