Abstract

   Toxocara canis, the causative agent of zoonotic toxocariasis in humans,
   is a parasitic roundworm of canids with a complex life cycle. While
   macrocyclic lactones (MLs) are successful at treating adult T. canis
   infections when used at FDA-approved doses in dogs, they fail to kill
   somatic third-stage larvae. In this study, we profiled the
   transcriptome of third-stage larvae derived from larvated eggs and
   treated in vitro with 10 μM of the MLs ivermectin and moxidectin. We
   analyzed transcriptional changes in comparison with untreated control
   larvae. In ivermectin-treated larvae, we identified 608 differentially
   expressed genes (DEGs), of which 453 were upregulated and 155 were
   downregulated. In moxidectin-treated larvae, we identified 1413 DEGs,
   of which 902 were upregulated and 511 were downregulated. Notably, many
   DEGs were involved in critical biological processes and pathways
   including transcriptional regulation, energy metabolism, body structure
   and function, physiological processes such as reproduction,
   excretory/secretory molecule production, host-parasite response
   mechanisms, and parasite elimination. We also assessed the expression
   of known ML targets and transporters, including glutamate-gated
   chloride channels (GluCls), and ATP-binding cassette (ABC)
   transporters, subfamily B, with a particular focus on P-glycoproteins
   (P-gps). We present gene names for previously uncharacterized T. canis
   GluCl and transporter genes using phylogenetic analysis of nematode
   orthologs to provide uniform gene nomenclature. Our study revealed that
   the expression of two GluCls and eight ABCB genes, particularly five
   P-gps were significantly altered in response to ML treatment. Compared
   to controls, Tca-glc-3, Tca-avr-14, Tca-haf-10, and Tca-Pgp-13.2 were
   downregulated in ivermectin-treated larvae, while Tca-abcb7,
   Tca-Pgp-11.2, and Tca-Pgp-2 were downregulated in moxidectin-treated
   larvae. Conversely, Tca-haf-9, Tca-Pgp-11.3, and Tca-Pgp-16.3 were
   upregulated in moxidectin-treated larvae. These findings suggest that
   MLs broadly impact transcriptional regulation in T. canis larvae.

   Keywords: Dogs, Moxidectin, Ivermectin, Toxocara canis, Toxocariasis,
   RNA-Seq, Transcriptome, Glutamate-gated chloride channels, ATP-binding
   cassette transporters, mRNA

Graphical abstract

   [28]Image 1
   [29]Open in a new tab

Highlights

     * •
       Macrocyclic lactones impact transcriptional responses in Toxocara
       canis third-stage larvae.
     * •
       608 (ivermectin) and 1413 (moxidectin) differentially expressed
       genes in ML-treated larvae.
     * •
       Two GluCls and eight ABCB genes were altered in response to MLs.

1. Introduction

   Toxocara canis is a parasitic roundworm that infects canids. Zoonotic
   T. canis larval infection in humans, known as toxocariasis, is a
   significant public health concern ([30]Hotez and Wilkins, 2009). In the
   U.S., it is estimated that 5 %–20 % of the population (∼17–66 million
   people) have been exposed to Toxocara spp., with notably higher
   exposure among individuals from lower socioeconomic backgrounds and
   those residing in the South and Northeastern parts of the country
   ([31]Congdon and Lloyd, 2011; [32]Berrett et al., 2017). On a global
   scale, the estimated number of people exposed is 19 % (1.52 billion
   people) ([33]Rostami et al., 2019). The risk of infection with T. canis
   is especially high in developing countries where poor sanitation and
   animal contact are common ([34]Holland, 2017).

   Humans become infected by accidentally consuming larvated eggs of T.
   canis ([35]Glickman and Shofer, 1987). Toxocariasis in humans has
   various clinical manifestations, ranging from infections with mild
   symptoms such as fever and abdominal pain to severe visceral and/or
   ocular larva migrans ([36]Despommier, 2003). To effectively mitigate
   the impact of toxocariasis on human health ([37]Chen et al., 2018), it
   is crucial to gain a deeper understanding of the biological processes
   and pathways involved in the transmission and infection of T. canis.

   The life cycle of T. canis involves the development of adult worms in
   the small intestine of infected dogs. Infection in dogs typically
   begins by the ingestion of larvated T. canis eggs found in contaminated
   soil or on the fur of infected animals ([38]Holland, 2017). Ingested
   eggs hatch in the duodenum, penetrate the intestinal mucosa, and
   migrate to the liver and lungs ([39]Schnieder et al., 2011). In dogs
   under three months of age, larvae are coughed up from the lungs and
   swallowed, completing their life cycle to adulthood in the small
   intestines. Adult female worms shed thousands of eggs, which pass out
   in the feces and contaminate the surrounding environment. These eggs
   can survive in the environment for long periods of time, creating an
   increased exposure risk for humans ([40]Azam et al., 2012).

   Impacts of T. canis infection in dogs can vary depending on the
   severity of the infection and the age of the animal ([41]Greve, 1971;
   [42]Oshima, 1976; [43]Dubey, 1978; [44]Fahrion et al., 2008). In severe
   infections, adult T. canis infection can lead to weight loss, poor
   growth, vomiting, and diarrhea ([45]Schnieder et al., 2011). Adult T.
   canis infections can be successfully treated using several anthelmintic
   drug classes including tetrahydropyrimidines (pyrantel pamoate, 5 mg/kg
   ([46]Clark et al., 1991)), benzimidazoles (fenbendazole, 50 mg/kg
   ([47]NADA, 1983), or febantel, 25 mg/kg ([48]Lloyd and Gemmell, 1992))
   and macrocyclic lactones (MLs), such as ivermectin (200 μg/kg
   ([49]Heredia Cardenas et al., 2017)), moxidectin (2.5 mg/kg ([50]NADA,
   2006)), and milbemycin oxime (0.5 mg/kg ([51]Bowman et al., 1988)). In
   older dogs, larvae undergo somatic migration, commonly arresting in the
   skeletal muscle, kidneys, and liver ([52]Manhardt and Stoye, 1981;
   [53]Schnieder et al., 2011). Somatic larvae can remain dormant for long
   periods and reactivate in dams during the third trimester of pregnancy
   ([54]Yutuc, 1949; [55]Webster, 1958; [56]Koutz et al., 1966;
   [57]Soulsby, 1983). Larvae are transmitted in utero from the infected
   dam to developing neonates ([58]Burke and Roberson, 1985). Thus,
   somatic larvae in a dog's muscles and viscera remain a persistent
   source of infection. Arrested somatic larvae of T. canis in dogs cannot
   be killed by anthelmintic treatments, including MLs, when used at
   FDA-approved labeled doses and therefore remain arrested in the
   tissues, demonstrating a perceived tolerance to the therapeutic effects
   of these compounds ([59]Overgaauw, 1997).

   The mechanism behind the tolerance of T. canis somatic larvae to MLs is
   not fully understood. MLs act on glutamate-gated chloride channels
   (GluCls) in parasitic nematodes ([60]Arena et al., 1995;
   [61]Wolstenholme and Rogers, 2005). GluCls are part of the cys-loop
   ligand-gated ion channels (LGICs) superfamily and play an important
   role in signal transduction in the invertebrate nervous system
   ([62]Wolstenholme, 2012). GluCls consist of five heterogeneous or
   homogeneous subunits ([63]Lamassiaude et al., 2022). Downregulation in
   the transcription of GluCls is associated with ivermectin resistance in
   the parasitic nematode Haemonchus contortus ([64]Williamson et al.,
   2011; [65]Tuersong et al., 2022). It is thought that alterations in
   GluCl subunit expression lead to different ion channel stoichiometry
   thereby altering ion channel pharmacology ([66]Whittaker et al., 2017).
   However, neither this phenomenon nor transcriptional changes in T.
   canis GluCls in response to ML treatment have been studied.

   Besides alterations in the GluCl target sites, there are other
   non-target-site mechanisms that may protect nematodes from MLs. In
   several nematodes, Permeability glycoproteins (P-gps), which belong to
   the ATP-binding cassette (ABC) transporter subfamily B (ABCB1) are
   involved with ML efflux, which is responsible for phenotypic resistance
   ([67]Blackhall et al., 1998; [68]Prichard and Roulet, 2007;
   [69]Bourguinat et al., 2008; [70]Dicker et al., 2011; [71]De Graef et
   al., 2013; [72]Janssen et al., 2013, [73]2015; [74]Mani et al., 2016;
   [75]Whittaker et al., 2017; [76]Figueiredo et al., 2018; [77]Gerhard et
   al., 2020). In T. canis, Tca-Pgp-11 has been demonstrated to transport
   ivermectin and be inhibited by known P-gp inhibitors ([78]Jesudoss
   Chelladurai et al., 2021). Additionally, thirteen P-gp genes are
   expressed in various larval stages of T. canis and have been studied
   using qPCR ([79]Jesudoss Chelladurai et al., 2023).

   To better understand the mechanisms underlying the effect of MLs on T.
   canis, a comprehensive analysis of its transcriptome is essential.
   Currently, there is limited data on the transcriptomes of T. canis,
   with one previous study comparing transcription in adult and larval
   isolates ([80]Zhu et al., 2015) and another comparing transcription in
   adult male and female isolates ([81]Zhou et al., 2017). Since the
   infectious stage of T. canis is the third-stage larvae (L3), there is a
   critical need to understand the transcriptional effects of MLs on this
   stage. Additionally, L3 larvae from eggs serve as a surrogate for
   somatic L3s, which are difficult to isolate from the tissues of
   infected dogs.

   In this study, we assessed the differential expression of genes in the
   comprehensive transcriptome of hatched T. canis third-stage larvae
   treated with the MLs ivermectin or moxidectin. We performed gene
   ontology analyses on the differentially expressed genes (DEGs).
   Additionally, we analyzed the transcription of annotated GluCl and ABCB
   genes of T. canis in response to ivermectin and moxidectin and further
   examined their expression using quantitative real-time polymerase chain
   reaction (qPCR).

2. Materials and methods

2.1. Parasites

   Adult T. canis male and female worms were opportunistically obtained
   from a naturally infected dog that was presented to the Kansas State
   University (KSU) Veterinary Diagnostic Laboratory for necropsy. Eggs
   were isolated from the uterus of adult gravid female worms and
   incubated at 25°C for three weeks to allow larval development to L3s.
   Hatching was induced as previously described ([82]Ponce-Macotela et
   al., 2011); larvae were isolated and suspended in Roswell Park Memorial
   Institute (RPMI) −1640 medium (Gibco by Thermo Fisher Scientific, Grand
   Island, NY). Hatched larvae (n = ∼500 - 1000) were subjected to three
   conditions: incubation with 10 μM ivermectin (IVM, Thermo Fisher
   Scientific, Ward Hill, MA), 10 μM moxidectin (MOX, Sigma-Aldrich, St.
   Louis, MO), or no treatment (controls, RPMI-1640 only), at 37°C for
   12 h and in biological triplicates. There was no observed loss in
   larval motility following drug incubation compared to controls. Samples
   were frozen at −80 °C until RNA extraction.

2.2. RNA isolation and RNA-seq

   Larvae were homogenized and RNA was isolated from each sample using the
   standard Trizol reagent extraction protocol with Phasemaker tubes
   (Invitrogen, Carlsbad, CA) as recommended by the manufacturer. Total
   RNA was quantified using a Qubit RNA Broad Range Assay Kit (Life
   Technologies/Invitrogen, Eugene, OR). RNA-seq was performed at the
   CEZID Molecular and Cellular Biology Core at KSU. Briefly, stranded
   mRNA was purified with oligo (dT) magnetic beads and libraries were
   prepared with an Illumina® Stranded mRNA Prep Ligation kit (Illumina,
   San Diego, CA) with a minimum of 100 ng of RNA. A brief graphical
   representation of these methods can be found in [83]Fig. 1. Paired-end
   sequencing (2x75 cycles) was performed on an Illumina NextSeq 550
   sequencer using a NextSeq 500/550 High Output Kit v2.5 (Illumina, San
   Diego, CA). Raw data from this project is available through the NCBI
   SRA database (Accession numbers SAMN38303355 - SAMN38303369;
   Bioproject: PRJNA1041894).

Fig. 1.

   [84]Fig. 1
   [85]Open in a new tab

   Brief graphical depiction of the experimental methods for the
   transcriptomic study of Toxocara canis third-stage larvae exposed to
   ivermectin, moxidectin, or RPMI-1640 (controls). This figure was
   created using [86]Biorender.com.

2.3. Bioinformatic analysis

   Raw data was obtained in FASTQ format and analyzed using the STAR
   pipeline on Galaxy servers ([87]Galaxy, 2022). Briefly, FASTQ files
   were trimmed using Trimmomatic (Version 0.38) ([88]Bolger et al.,
   2014). Trimmed reads were aligned and mapped to the genes annotated in
   the T. canis draft genome (NCBI GenBank assembly GCA_000803305.1), and
   gene expression was measured as reads per gene using STAR ([89]Dobin et
   al., 2013). Raw STAR data logs can be found in [90]Supplemental
   Document 1. Raw and trimmed read qualities were assessed with FastQC
   (Version 0.12.1). Alignment quality was assessed with Qualimap BAMQC
   (Version 2.2.2c) ([91]Garcia-Alcalde et al., 2012; [92]Okonechnikov et
   al., 2016). Reads per gene counts from STAR were analyzed using custom
   R scripts (Version 4.3.2). Gene expression normalization, and
   differential expression analysis in response to drug treatment were
   performed using DESeq2 (Version 1.40.1) ([93]Love et al., 2014). DESeq2
   performs negative binomial Wald test by Benjamini–Hochberg correction
   for multiple comparisons. Differentially expressed genes (DEGs) were
   defined as those with a false discovery rate (FDR)–adjusted p-value <
   0.05 and an absolute log2 fold change ≥ 1 for upregulated genes or
   log2fold change < −1 for downregulated genes. Principal components
   analysis (PCA) was done on variance-stabilizing transformed (vst)
   normalized values. Volcano plots were created with EnhancedVolcano
   (Version 1.18.0) ([94]Blighe et al., 2018). Pearson correlations as a
   measure of intragroup invariance were calculated and DEGs were filtered
   with tidyverse (Version 2.0.0). Venn diagrams to determine DEGs that
   were shared between drug treatments were created with ggvenn (Version
   0.1.10). Functional gene ontology analysis was plotted with gprofiler2
   (Version 0.2.2) ([95]Kolberg et al., 2020). Gene set enrichment on DEGs
   was performed in ShinyGO (Version 0.80) ([96]Ge et al., 2020).

   Gene-specific information about 25 ML-associated genes (transporters
   and targets) was obtained using rentrez (Version 1.2.3) ([97]Winter,
   2017). Phylogenetic analysis was performed using NGPhylogeny.fr
   ([98]Lemoine et al., 2018, [99]2019). Briefly, protein sequences were
   manually curated from the Genbank protein database, aligned with MAFFT
   ([100]Katoh and Standley, 2013), and trimmed with BMGE ([101]Criscuolo
   and Gribaldo, 2010). A maximum likelihood (ML) tree was constructed
   with PhyML 3.0 ([102]Guindon et al., 2010) after determining the best
   model using SMS ([103]Lefort et al., 2017). The tree in Newick format
   was visualized using iTOL (version 6) ([104]Junier and Zdobnov, 2010;
   [105]Letunic and Bork, 2024) and ggtree (Version 3.10.1) ([106]Xu et
   al., 2022). Heatmaps were created with variance stabilized and
   transformed DESeq2 data using ComplexHeatmap (Version 2.16.0) ([107]Gu
   et al., 2016), ggpubr (Version 0.6.0), and grid (Version 4.3.2).
   Boxplots to visualize expression levels were created using ggplot2
   (Version 3.5.0). Differences among treatments were first assessed using
   the Kruskal–Wallis test, followed by Dunn's post-hoc test with rstatix
   (Version 0.7.2). Other R libraries used included ggprism (Version
   1.0.5) and patchwork (Version 1.2.0).

2.4. qPCR analysis of larvae following in vitro drug exposure

   Expression levels of a subset of ML-associated genes (15 genes) were
   determined using qPCR. qPCR reactions were individually optimized using
   diluted cDNA synthesized from adult T. canis worms. Specificity was
   determined using melt curve analysis and sequencing. Primers used in
   this study can be found in [108]Supplemental Table 1. 18S rRNA was used
   as a reference gene and amplified using previously described primers
   ([109]Jesudoss Chelladurai et al., 2023). Larval cDNA was synthesized
   from 100 ng of total RNA in a 20 μL volume using the High-Capacity cDNA
   Reverse Transcription Kit (Applied Biosystems) and stored at −20 °C
   till use. To determine differences between drug exposed and unexposed
   larvae, qPCR was carried out in biological triplicates and technical
   duplicates in a volume of 20 μL with 2 μL of diluted cDNA, PowerTrack
   SYBR Green Master Mix (Applied Biosystems), and 200–500 nM of diluted
   primers. PCR efficiency was determined for each primer pair using
   LinRegPCR (Version 2021.2) ([110]Ramakers et al., 2003). For qPCR
   analyses, relative expression was calculated using the Pfaffl method
   ([111]Pfaffl and Bustin, 2004). Boxplots to visualize expression levels
   were created using ggplot2 (Version 3.5.0). Differences among
   treatments were assessed by the Kruskal–Wallis test, followed by Dunn's
   post-hoc test with rstatix (Version 0.7.2).

3. Results

3.1. Sequencing

   We obtained transcriptomic sequences from hatched L3s that were
   untreated, or treated with ivermectin, or treated with moxidectin
   (n = 3, ∼500–1000 larvae each) using an Illumina NextSeq 550 sequencer.
   A total of 295,503,166 paired-end raw reads were generated, with an
   average of 32,833,685 reads per sample. After adapter trimming, a total
   of 293,127,407 pair-end reads were retained, with an average of
   32,569,712 reads per sample. The reads were assembled against the
   annotated reference genome (GenBank Accession #PRJNA248777; 18,596
   genes ([112]Zhu et al., 2015)) before transcripts were quantified using
   STAR. The sequencing data generated was of high quality, with an
   average of 81.3 % of the reads mapping to the T. canis draft genome
   ([113]Table 1).

Table 1.

   Summary of the features of the sequenced cDNA libraries for samples
   from this study.
   Sample name Sample group Paired read count (raw) Nucleotide count (raw
   Paired read count (trimmed) Nucleotide count (trimmed) Mapped read
   count (%)
   Control T1 Control T. canis L3 32,249,245 2.3 Gbp 32,174,441 2.2 Gbp
   84.54 %
   Control T7 Control T. canis L3 54,071,313 3.9 Gbp 54,010,559 3.8 Gbp
   86.05 %
   Control T16 Control T. canis L3 40,128,609 2.9 Gbp 40,087,855 2.8 Gbp
   87.9 %
   Ivermectin T3 Ivermectin T. canis L3 22,121,380 1.6 Gbp 21,543,247 1.4
   Gbp 81.79 %
   Ivermectin T4 Ivermectin T. canis L3 17,179,802 1.2 Gbp 17,105,795 1.2
   Gbp 92.92 %
   Ivermectin T16 Ivermectin T. canis L3 26,457,311 1.9 Gbp 25,746,262 1.8
   Gbp 91.54 %
   Moxidectin T2 Moxidectin T. canis L3 29,364,220 2.1 Gbp 29,186,760 2.0
   Gbp 92.84 %
   Moxidectin T7 Moxidectin T. canis L3 36,440,579 2.6 Gbp 36,389,622 2.5
   Gbp 45.54 %
   Moxidectin T8 Moxidectin T. canis L3 37,490,707 2.7 Gbp 36,882,866 2.6
   Gbp 68.59 %
   [114]Open in a new tab

   To determine similarities between the untreated and treated larval
   groups a PCA was conducted to visualize the relationship between the
   groups based on vst of normalized expression values ([115]Fig. 2A).
   Principal component 1 explained 39 % of the variance observed.
   Principal component 2 explained 29 % of the variance observed.
   Intragroup Pearson correlation within the group ranged from 0.94 to
   0.99 for the control group, 0.71 to 0.95 for the ivermectin-treated
   group, and 0.97 to 0.98 for the moxidectin-treated group. Untreated
   controls clustered separately from the ivermectin- and
   moxidectin-treated larvae, indicating that the drug treatments had an
   impact on the gene expression profiles of the larvae.

Fig. 2.

   [116]Fig. 2
   [117]Open in a new tab

   (A) Principal Component Analysis (PCA) of Toxocara canis third-stage
   larvae (L3s) exposed to RPMI-1640 medium (controls; blue circles),
   ivermectin (yellow squares), or moxidectin (red diamonds). Each point
   represents a pool of ∼500 - 1000 L3s (one biological replicate); (B–C)
   Volcano plots of differentially expressed genes in (B) ivermectin- or
   (C) moxidectin-treated larvae compared to untreated controls.
   Significance thresholds were set at log2FC > 1 (two-fold change) for
   upregulated genes, and log2FC < −1 for downregulated genes and adjusted
   p-value < 0.05. Upregulated genes are shown in red, downregulated genes
   in blue, and non-significant genes in gray. A total of 18,596 genes
   were analyzed, with 454 significantly upregulated and 155 significantly
   downregulated in ivermectin-treated L3s and 902 significantly
   upregulated and 511 significantly downregulated in moxidectin-treated
   L3s; (D–E) Venn diagrams depicting the overlap of shared upregulated
   (D) and downregulated (E) genes in larvae treated with ivermectin and
   moxidectin. The upregulated genes show an overlap of 240 shared genes,
   while the downregulated genes show an overlap of 83 shared genes; (F)
   Heatmap of normalized expression (vst) values for all 18,596 genes
   annotated in the T. canis genome for control, ivermectin-, and
   moxidectin-treated L3s. Color intensity represents normalized
   expression values, with the scale ranging from black (lowest
   expression) to red (highest expression).

3.2. Analysis of differentially expressed transcripts

   To identify differentially expressed transcripts, we used the
   reference-guided transcript assembly method using STAR and DESeq2.
   Volcano plots showing fold changes and adjusted p-values for the 18,596
   genes are shown in [118]Fig. 2B–C. We identified 608 DEGs in larvae
   treated with ivermectin compared to controls, of which 453 were
   upregulated (two-fold change, that is, log2fold change > 1 and adjusted
   p-value < 0.05), and 155 were downregulated (log2fold change < −1 and
   adjusted p-value < 0.05), respectively ([119]Fig. 2B). Additionally, we
   identified 1413 DEGs in larvae treated with moxidectin compared to
   controls, of which 902 were upregulated (log2fold change > 1 and
   adjusted p-value < 0.05), and 511 were downregulated (log2fold
   change < −1 and adjusted p-value < 0.05) ([120]Fig. 2C). We identified
   240 upregulated genes ([121]Figs. 2D), and 83 downregulated genes
   ([122]Fig. 2E) shared between the ivermectin and moxidectin-treated
   larvae.

   To visualize the expression patterns of the differentially expressed
   transcripts in ivermectin- and moxidectin-treated larvae compared to
   controls, a heatmap of all expressed genes was generated and was sorted
   by gene expression ([123]Fig. 2F). We observed that approximately
   11.22 % of genes exhibited high expression levels, with normalized
   values between 10 and 15, indicated in yellow (10.94 %), and between 15
   and 25, indicated in red (0.28 %). Conversely, a significant number of
   genes showed little to no expression, indicated in black (30 %). The
   hierarchical clustering revealed two predominate clusters across
   different conditions, with further hierarchical sub-clustering within
   each.

   To gain insight into the pathways and functions significantly
   associated with the observed changes in gene expression, a gene
   ontology (GO) enrichment analysis was performed on the DEGs ([124]Fig.
   3, [125]Fig. 4, [126]Fig. 5). DEGs were associated with GO terms
   categorized as biological processes (BP), molecular functions (MF), and
   cellular components (CC) and described in the following sections.

Fig. 3.

   [127]Fig. 3
   [128]Open in a new tab

   (A&C) Gene Ontology (GO) bubble plots of genes (A) upregulated and (C)
   downregulated in Toxocara canis third-stage larvae (L3s) treated with
   ivermectin. The size of the bubble represents the number of genes
   associated with each GO term, and the color intensity of the bubble
   indicates the significance (adjusted p-value) of enrichment. GO terms
   are grouped by biological process (BP, blue), molecular function (MF,
   red), and cellular component (CC, yellow) categories; (B&D) Lollipop
   plots representing pathway enrichment analysis for the top 10 most
   significantly enriched pathways for genes (B) upregulated and (C)
   downregulated in T. canis ivermectin-treated L3s. The length and color
   of the lollipops correspond to fold enrichment, with blue indicating
   lower fold enrichment and red indicating higher fold enrichment. The
   size of the dots represents the number of genes associated with each
   pathway.

Fig. 4.

   [129]Fig. 4
   [130]Open in a new tab

   (A&C) Gene Ontology (GO) bubble plots of genes (A) upregulated and (C)
   downregulated in Toxocara canis third-stage larvae (L3s) treated with
   moxidectin. The size of the bubble represents the number of genes
   associated with each GO term, and the color intensity of the bubble
   indicates the significance (adjusted p-value) of enrichment. GO terms
   are grouped by biological process (BP, blue), molecular function (MF,
   red), and cellular component (CC, yellow) categories; (B&D) Lollipop
   plots representing pathway enrichment analysis for the top 10 most
   significantly enriched pathways for genes (B) upregulated and (C)
   downregulated in T. canis moxidectin-treated L3s. The length and color
   of the lollipops correspond to fold enrichment, with blue indicating
   lower fold enrichment and red indicating higher fold enrichment. The
   size of the dots represents the number of genes associated with each
   pathway.

Fig. 5.

   [131]Fig. 5
   [132]Open in a new tab

   (A&C) Gene Ontology (GO) bubble plots of the overlapping (A)
   upregulated and (C) downregulated shared genes between ivermectin- and
   moxidectin-treated Toxocara canis third-stage larvae (L3s). The size of
   the bubble represents the number of genes associated with each GO term,
   and the color intensity of the bubble indicates the significance
   (adjusted p-value) of enrichment. GO terms are grouped by biological
   process (BP, blue), molecular function (MF, red), and cellular
   component (CC, yellow) categories; (B&D) Lollipop plots representing
   pathway enrichment analysis for the top 10 most significantly enriched
   pathways for the overlapping (B) upregulated and (C) downregulated
   shared genes between ivermectin- and moxidectin-treated T. canis L3s.
   The length and color of the lollipops correspond to fold enrichment,
   with blue indicating lower fold enrichment and red indicating higher
   fold enrichment. The size of the dots represents the number of genes
   associated with each pathway.

3.3. Response to ivermectin

   The 453 T. canis genes upregulated in response to ivermectin exposure
   are listed in [133]Supplemental Table 2A. Among the most upregulated
   genes in response to ivermectin exposure were a mediator of RNA
   polymerase II transcription subunit 13, Tca-let-19 (Tcan_05892,
   log2fold change: 5, padj. 5.10E-03), glutamate dehydrogenase 2,
   mitochondrial, Tca-GLUD2 (Tcan_03811, log2fold change: 3.9, padj.
   2.32E-02), nuclear anchorage protein 1, Tca-anc-1 (Tcan_00659, log2fold
   change: 3.5, padj. 1.09E-04), putative U5 small nuclear
   ribonucleoprotein helicase, Tca-Y46G5.4 (Tcan_05430, log2fold change:
   3.4, padj. 3.43E-03); additionally, there were six hypothetical
   proteins (Tcan_00335, Tcan_06163, Tcan_07262, Tcan_15031, Tcan_00594,
   and Tcan_00011), of which Tcan_06163 was identified as an ortholog of
   C. elegans T20G5.12.

   Gene ontology enrichment analysis of the upregulated genes for larvae
   exposed to ivermectin ([134]Supplemental Table 3A/[135]Fig. 3A and B)
   revealed a significant enrichment of genes with GO terms associated
   with “polymorphism” (FDR = 6.6E-04), “nucleosome positioning”
   (FDR = 6.6E-04), “mixed, incl. troponin i binding, and striated muscle
   myosin thick filament” (FDR = 1.7E-04), “mostly uncharacterized, incl.
   serine-type endopeptidase inhibitor activity” (FDR = 4.3E-05), “mixed,
   incl. troponin complex, and sarcoplasmic reticulum” (FDR = 1.3E-08),
   ”respiratory chain” (FDR = 3.7E-04), “mostly uncharacterized, incl.
   troponin i binding, and rna destabilization” (FDR = 46.6E-04),
   "striated muscle thin filament” (FDR = 42.7E-05) “myofilament”
   (FDR = 2.2E-06), and “NADH dehydrogenase (ubiquinone) activity”
   (FDR = 6.6E-04).

   The 155 T. canis genes downregulated in response to ivermectin exposure
   are listed in [136]Supplemental Table 2B. Among the most downregulated
   genes in response to ivermectin were a chondroitin proteoglycan 4,
   Tca-cpg-4 (Tcan_10605, log2fold change: 4, padj. 3.76E-03),
   methylmalonyl-CoA epimerase, mitochondrial, Tca-Mcee (Tcan_12863,
   log2fold change: 3.3, padj. 2.18E-02), metalloproteinase inhibitor 2,
   Tca-TIMP2 (Tcan_03888, log2fold change: 3.3, padj. 3.28E-02), putative
   oxidoreductase, Tca-dhs-27 (Tcan_13635, log2fold change: 3.3, padj.
   3.88E-02), amidophosphoribosyltransferase, Tca-Prat (Tcan_18524,
   log2fold change: 3.0, padj. 2.32E-03), nanos -like protein 2,
   Tca-Nanos2 (Tcan_18885, log2fold change: 3.0, padj. 1.21E-02), and
   UDP-glucuronosyltransferase 2C1, Tca-UGT2C1 (Tcan_05445, log2fold
   change: 2.9, padj. 6.19E-06). There were three hypothetical proteins
   (Tcan_15481, Tcan_14498, and Tcan_01557).

   Gene ontology enrichment analysis of the downregulated genes for larvae
   exposed to ivermectin ([137]Supplemental Table 3B/[138]Fig. 3C and D)
   revealed a significant enrichment of genes with GO terms associated
   with “npBAF complex” (FDR = 2.1E-02), “amidophosphoribosyltransferase
   activity” (FDR = 2.8E-02), “protein retention in er lumen”
   (FDR = 2.8E-02), “phosphoribosylaminoimidazole carboxylase complex”
   (FDR = 2.8E-02), “mixed, incl. ultradian rhythm, and kelch repeat”
   (FDR = 3.5E-02), “amidophosphoribosyltransferase activity, and
   l-glutamine aminotransferase activity” (FDR = 3.5E-02), “mixed, incl.
   gamma-tubulin binding, and fatty-acyl-coa binding” (FDR = 3.5E-02),
   “maintenance of protein localization in endoplasmic reticulum”
   (FDR = 3.5E-02), “lysosome, and transmembrane transporter activity”
   (FDR = 4.6E-02), and “protein folding chaperone, and heat shock protein
   binding” (FDR = 4.6E-02).

3.4. Response to moxidectin

   The 902 T. canis genes upregulated in response to moxidectin exposure
   are listed in [139]Supplemental Table 4A. Among the most upregulated
   genes in response to moxidectin were a trypsin inhibitor, partial
   (Tcan_08038, log2fold change: 5.6, padj. 1.16E-04), glutamate
   dehydrogenase 2, mitochondrial, Tca-GLUD2 (Tcan_03811, log2fold change:
   4.9, padj. 1.12E-03), group XV phospholipase A2, Tca- PLA2G15
   (Tcan_15544, log2fold change: 4.7, padj. 7.56E-03), aspartyl protease
   inhibitor, partial, Tca-API-3 (Tcan_02586, log2fold change: 4.3, padj.
   9.96E-04), and protocadherin-16, Tca-F32A5.4 (Tcan_07881, log2fold
   change: 3.8, padj. 3.56E-09). Five hypothetical proteins were also
   found (Tcan_01052, Tcan_00040, Tcan_15628, Tcan_00055, Tcan_00043), of
   which Tcan_00043 was identified as an ortholog to V-type proton ATPase
   subunit a in Dirofilaria immitis or Bm33 in Brugia malayi ([140]Ghedin
   et al., 2007; [141]Gomes-de-Sa et al., 2022).

   Gene ontology enrichment analysis of the upregulated genes for larvae
   exposed to moxidectin ([142]Supplemental Table 5A/[143]Fig. 4A and B)
   revealed a significant enrichment of genes with GO terms associated
   with “ribosomal protein” (FDR = 1.31E-38), “cytosolic large ribosomal
   subunit” (FDR = 9.45E-34), “cytosolic large ribosomal subunit, and
   cytosolic small ribosomal subunit” (FDR = 3.43E-45), “cytosolic
   ribosome, and preribosome, small subunit precursor” (FDR = 5.83E-53),
   “cytosolic ribosome” (FDR = 2.01E-48), “cytosolic small ribosomal
   subunit” (FDR = 1.45E-16), “cytosolic ribosome, and eukaryotic
   translation initiation factor 3 complex” (FDR = 3.65E-51), “NADH
   dehydrogenase activity” (FDR = 7.15E-13), “ribonucleoprotein”
   (FDR = 9.65E-34), and “cytosolic ribosome, and translational
   initiation” (FDR = 1.58E-49).

   The 511 T. canis genes downregulated in response to moxidectin exposure
   are listed in [144]Supplemental Table 4B. Among the most downregulated
   genes in response to moxidectin were an arginase, Tca-rocF (Tcan_06332,
   log2fold change: 4.2. padj. 2.42E-02), cuticlin-1, Tca-cut-1
   (Tcan_11029, log2fold change: 3.5, padj. 1.28E-02), alpha-aminoadipic
   semialdehyde synthase, mitochondrial, Tca-AASS (Tcan_18345, log2fold
   change: 3.4, padj. 1.32E-06), TRIO and F-actin-binding protein,
   partial, Tca- Triobp (Tcan_10839, log2fold change: 2.8, padj.
   9.43E-03), lysine histidine transporter 1, Tca-LHT1 (Tcan_17795,
   log2fold change: 2.8, padj. 1.20E-03), putative zinc finger protein
   F56D1.1, partial, Tca-F56D1.1 (Tcan_16300, log2fold change: 2.6, padj.
   3.47E-03), and methylmalonic aciduria and homocystinuria type D -like
   protein, mitochondrial, Tca- Mmadhc (Tcan_06175, log2foldchange: 2.5,
   padj.1.43E-04). Three hypothetical proteins were identified
   (Tcan_04808, Tcan_02388, and Tcan_09119).

   Gene ontology enrichment analysis of the downregulated genes for larvae
   exposed to moxidectin ([145]Supplemental Table 5B/[146]Fig. 4C and D)
   revealed a significant enrichment of genes with GO terms associated
   with “phosphatidylethanolamine biosynthetic process”
   (FDR = 0.00037204), “phosphatidylethanolamine metabolic process”
   (FDR = 0.000824851), “endoplasmic reticulum-golgi intermediate
   compartment” (FDR = 0.001363654), “recycling endosome”
   (FDR = 0.000785154), “mixed, incl. cellular response to hormone
   stimulus, and longevity regulating pat” (FDR = 0.001076538), “molecular
   adaptor activity” (FDR = 0.000374112), “mixed, incl. phosphatase
   complex, and hippo signaling pathway – multiple species”
   (FDR = 0.001056856), “golgi membrane” (FDR = 0.001403698), “cytoplasmic
   vesicle” (FDR = 6.59E-08), and “intracellular vesicle”
   (FDR = 6.59E-08).

3.5. Shared analysis

   To identify DEGs shared between ivermectin- and moxidectin-treated
   larvae, we performed a comparative analysis. DEGs shared in response to
   ivermectin and moxidectin exposure are listed in [147]Supplemental
   Table 6. Among the top shared upregulated genes were glutamate
   dehydrogenase 2, mitochondrial, Tca-GLUD2 (Tcan_03811, log2fold change:
   4.4), protocadherin-16, Tca-F32A5.4 (Tcan_07881, log2fold change: 3.0),
   cuticle collagen 40, Tca-col-40 (Tcan_09436, log2fold change: 2.8), and
   two hypothetical proteins (Tcan_14010 and Tcan_15628). Tcan_14010 was
   identified as an ortholog to a “neo-calmodulin-like” protein in
   Ylistrum balloti (Ballot's saucer scallop). Gene ontology enrichment
   analysis of upregulated shared genes ([148]Supplemental Table
   7A/[149]Fig. 5A and B) revealed a significant enrichment of genes with
   GO terms associated with “troponin complex” (FDR = 3.7E-04),
   “respiratory chain” (FDR = 3.7E-04), “mitochondrial electron transport,
   nadh to ubiquinone” (FDR = 9.3E-05), “NADH dehydrogenase (ubiquinone)
   activity” (FDR = 2.7E-05), and “Oxidoreductase activity, acting on
   nad(p)h, quinone or similar compound as acceptor” (FDR = 3.5E-05).

   Among the top shared downregulated genes were mitochondrial
   sodium/hydrogen exchanger 9B2, Tca-SLC9B2 (Tcan_12528, log2fold change:
   2.5), methylmalonic aciduria and homocystinuria type D -like protein,
   mitochondrial, Tca-Mmadhc (Tcan_06175, log2fold change: 2.3), T-cell
   activation inhibitor, mitochondrial, Tca-TCAIM (Tcan_14999, log2fold
   change: 2.3), as well as two hypothetical proteins (Tcan_00180 and
   Tcan_01557). Gene ontology enrichment analysis of downregulated shared
   genes ([150]Supplemental Table 7B/[151]Fig. 5C and D) revealed a
   significant enrichment of genes with GO terms associated with “npBAF
   complex” (FDR = 8.1E-03), “lysosome, and transmembrane transporter
   activity” (FDR = 1.8E-02), “protein folding chaperone, and heat shock
   protein binding” (FDR = 1.8E-02), “polypeptide
   n-acetylgalactosaminyltransferase activity” (FDR = 5.0E-03), and
   “nucleosome disassembly” (FDR = 2.3E-02).

3.6. Genes associated with ML function and transport

   To understand transcriptional changes in response to ML exposure in T.
   canis larvae, we identified 25 genes of interest annotated in the
   genome as ABCB transporters or GluCls and obtained their normalized
   expression (vst) values from the dataset. ABCB1 transporters of T.
   canis had been named previously ([152]Jesudoss Chelladurai et al.,
   2023) based on phylogenetic analysis following the recommendation of
   ([153]Beech et al., 2010).

   GluCls and ABCB transporters other than P-gps of T. canis have not been
   studied, and gene nomenclature has not been adequately addressed. To
   clarify gene naming for T. canis GluCls, we identified seven GluCls in
   the transcriptome of T. canis and assigned gene names to them
   ([154]Table 2) using a maximum likelihood phylogenetic tree ([155]Fig.
   6A). The phylogenetic analysis showed that orthologs of GluCls in T.
   canis clustered closely with those of other described nematodes,
   forming well-supported clades. Additionally, we identified one T. canis
   gene (GenBank Accession [156]KHN83224) that was annotated as a GluCl;
   however, phylogenetic analysis revealed that it clustered within a
   distinct clade alongside other nematode ligand-gated ion channels.
   Similarly, we annotated the other ABCB transporter family members,
   which were composed of half ATP-binding cassette transporters (haf)
   similar to those described in C. elegans ([157]Fig. 6B–[158]Table 3).
   Thus, we included 18 ABCB transporter genes, and seven GluCls in our
   analysis.

Table 2.

   Nomenclature for GluCl protein sequences in Toxocara canis.
   GenBank accession number (Protein database) Assigned gene names
   [159]KHN83814                               Tca-glc-2
   [160]KHN87150                               Tca-glc-5.2
   [161]KHN83224                               Tca-lgc
   [162]KHN82095                               Tca-glc-4.1
   [163]KHN81457                               Tca-avr-14
   [164]KHN80184                               Tca-glc-3
   [165]KHN77193                               Tcaglc-5.1
   [166]KHN74877                               Tca-glc-4.2
   [167]Open in a new tab

Fig. 6.

   [168]Fig. 6
   [169]Open in a new tab

   A maximum-likelihood phylogenetic tree of Toxocara canis
   glutamate-gated chloride channels (GluCls) (A) and half ATP-binding
   cassette (haf) transporters (B) with related nematode sequences
   obtained from GenBank. Both trees were constructed using PhyML with the
   SMS to determine the best fitting model and visualized using iTOL and
   ggtree. Branch support values were calculated from 1000 bootstrap
   replicates, with values > 75 shown. The GluCl tree (A) is rooted using
   a GluCl sequence from Aplysia californica (California sea hare) and the
   haf tree (B) was rooted using a P-gp sequence from Canis lupus
   familiaris (dog).

Table 3.

   Nomenclature for ABCB protein sequences in Toxocara canis.
   GenBank accession number (Protein database) Assigned gene names Full or
   partial Reference
   [170]KHN80157 Tca-Pgp-2 Full [171]Jesudoss Chelladurai et al. (2023)
   [172]KHN78383 Tca-Pgp-3 Partial [173]Jesudoss Chelladurai et al. (2023)
   [174]KHN77280 Tca-Pgp-9.1 Partial [175]Jesudoss Chelladurai et al.
   (2023)
   [176]KHN76604 Tca-Pgp-9.2 Partial [177]Jesudoss Chelladurai et al.
   (2023)
   [178]KHN84692 Tca-Pgp-10 Full [179]Jesudoss Chelladurai et al. (2023)
   [180]KHN73709 Tca-Pgp-11.1 Full [181]Jesudoss Chelladurai et al. (2023)
   [182]KHN87227 Tca-Pgp-11.2 Full [183]Jesudoss Chelladurai et al. (2023)
   [184]KHN89031 Tca-Pgp-11.3 Full [185]Jesudoss Chelladurai et al. (2023)
   [186]KHN87070 Tca-Pgp-13.1 Full [187]Jesudoss Chelladurai et al. (2023)
   [188]KHN87068 Tca-Pgp-13.2 Full [189]Jesudoss Chelladurai et al. (2023)
   [190]KHN80315 Tca-Pgp-16.1 Full [191]Jesudoss Chelladurai et al. (2023)
   [192]KHN80341 Tca-Pgp-16.2 Full [193]Jesudoss Chelladurai et al. (2023)
   [194]KHN86334 Tca-Pgp-16.3/3.2 Full [195]Jesudoss Chelladurai et al.
   (2023)
   [196]KHN88383 Tca-haf-6 Partial Present study
   [197]KHN83466 Tca-abcb7 Partial Present study
   [198]KHN73719 Tca-haf-10 Partial Present study
   [199]KHN73722 Tca-haf-4 Partial Present study
   [200]KHN76663 Tca-haf-9 Partial Present study
   [201]Open in a new tab

   We analyzed the expression of the identified GluCls and ABCBs
   ([202]Jesudoss Chelladurai et al., 2023) ([203]Fig. 7A). In our initial
   differential expression analysis of all 18,596 genes, GluCls were not
   found among the DEGs and only Tca-haf-9 was upregulated in the
   moxidectin-treated larvae compared to controls (Tcan_14960, log2fold
   change: 1.314). Upon further analysis using Krusal-Wallis tests with
   Dunn's multiple comparisons tests, we observed differences in the
   expression of nine out of the 25 genes of interest ([204]Fig. 7B). The
   Dunn's’s test identified that Tca-glc-3 was significantly downregulated
   in ivermectin-treated larvae compared to controls (p = 0.025). Of the
   ABCBs, compared to controls, Tca-abcb7, Tca-Pgp-11.2, and Tca-Pgp-2
   were downregulated in moxidectin-treated larvae (p = 0.025, 0.025, and
   0.037, respectively). Tca-haf-10 and Tca-Pgp-13.2 were downregulated in
   ivermectin-treated larvae compared to controls (p = 0.025 and 0.025,
   respectively). Conversely, Tca-Pgp-11.3, Tca-Pgp-16.3, Tca-haf-9 were
   upregulated in moxidectin-treated larvae compared to controls
   (p = 0.011, p = 0.025, and p = 0.017, respectively).

Fig. 7.

   [205]Fig. 7
   [206]Open in a new tab

   (A) Heatmap of normalized expression (vst) values for all 25 Toxocara
   canis genes of interest, including seven glutamate-gated chloride
   channels (GluCls), and 19 ATP-Binding Cassette Transporters, subfamily
   B (13 of which are P-glycoproteins (Pgps)), for control, ivermectin-,
   and moxidectin-treated third-stage larvae (L3s). Color intensity
   represents normalized expression values, with the scale ranging from
   black (lowest expression) to red (highest expression). (B) Boxplots of
   normalized expression (vst) values for the 25 T. canis genes of
   interest. The y-axis represents normalized expression. Statistical
   significance between treatment groups was assessed using the
   Kruskal-Wallis test followed by Dunn's post-hoc comparisons. Control
   L3s are represented as blue, ivermectin-treated L3s as yellow, and
   moxidectin-treated L3s as red shaded boxes.

   To confirm our transcriptomic findings, we performed qPCR on 15 genes
   of interest, while the remaining genes had been previously examined
   under ML treatment ([207]Jesudoss Chelladurai et al., 2023). This
   analysis showed that Tca-avr-14 and Tca-Pgp-13.2 were significantly
   downregulated in ivermectin-treated larvae compared to controls
   (p = 0.025, and p = 0.017, respectively; [208]Fig. 8). Additionally,
   compared to moxidectin-treated larvae, Tca-glc-2 was downregulated in
   ivermectin-treated larvae (p = 0.02).

Fig. 8.

   [209]Fig. 8
   [210]Open in a new tab

   qPCR validation of RNA-Seq results. Fold change in expression of
   15 ML-associated genes of interest following ivermectin (yellow) or
   moxidectin (red) treatment of larvae hatched in vitro compared to
   control larvae (blue). Fold change was calculated using the Pfaffl
   method using T. canis 18S rRNA as the reference gene. Statistical
   significance was assessed using the Kruskal–Wallis test followed by
   Dunn's post-hoc comparisons.

4. Discussion

   T. canis infections remain a persistent threat to both humans and
   animals globally, despite the availability of effective pharmaceuticals
   for treating adult infections. Somatic larvae within canine hosts serve
   as the primary reservoir for propagating the life cycle and are not
   cleared using MLs at approved doses ([211]Overgaauw, 1997). Similarly,
   hatched third-stage larvae observed in this study were also refractory
   to ML-induced killing, highlighting a shared tolerance across
   developmental stages. The mechanisms behind the tolerance of T. canis
   larvae to MLs are still poorly understood. In this study, we assessed
   the effects of the MLs on T. canis third-stage larvae in vitro using
   transcriptomics, focusing on DEGs and the expression levels of genes
   associated with ML action and resistance in nematodes, specifically
   GluCls and P-gps.

   To increase the rigor of the study, we used three replicates of larvae
   hatched from eggs from different adult worms and subjected them to
   10 μM ivermectin, moxidectin, or no drugs (RPMI-1640) for 12 h and
   performed RNA-seq on the Illumina platform. We observed a high
   intragroup Pearson correlation of mRNA transcripts in the control group
   (0.94–0.99) indicating strong homogeneity among the samples. The
   ivermectin group showed greater variability (0.71–0.95), while the
   moxidectin group (0.97–0.98) demonstrated a consistency similar to the
   control group. Furthermore, we observed significant differential
   expression of genes in larvae exposed to ivermectin (609 genes) and
   moxidectin (1413 genes) compared to controls. Of these, 453 were
   upregulated and 155 were downregulated in response to ivermectin
   exposure, while 902 were upregulated and 511 were downregulated in
   response to moxidectin exposure. Among these, 240 were shared between
   the two MLs in the upregulated genes and 83 in the downregulated gene
   set.

   The distinct transcriptional responses we observed in T. canis larvae
   following ivermectin and moxidectin exposure may reflect underlying
   differences in their molecular interactions. Although both belong to
   the macrocyclic lactone class, ivermectin is an avermectin, whereas
   moxidectin is a milbemycin. Their structural differences have been
   shown to influence pharmacological properties including lipophilicity,
   target affinity, and substrate specificity for P-gps ([212]Prichard et
   al., 2012). In C. elegans, ivermectin and moxidectin produced distinct
   phenotypic and transcriptional responses, and showed differential
   reliance on specific GluCl subunits ([213]Ardelli et al., 2009). In
   addition, Gerhard et al. showed that ivermectin uptake is strongly
   enhanced by pharyngeal pumping, while moxidectin activity is not,
   suggesting that moxidectin uptake may rely more heavily on absorption
   through the cuticle ([214]Gerhard et al., 2021). These differences may
   ultimately drive transcriptional changes in nematodes in response to
   exposure to these MLs.

   The consensus site of action of the MLs are the GluCls. GluCls are
   pentameric channels composed of heterogeneous or homogeneous subunits.
   MLs bind to these channels between the M1 and M3 transmembrane domains,
   specifically at the Leu217 and Ile222 sites on M1 ([215]Hibbs and
   Gouaux, 2011). GluCls have been extensively studied in Ascaris suum
   ([216]Jagannathan et al., 1999), C. elegans ([217]Arena et al., 1992;
   [218]Cully et al., 1994, [219]1996; [220]Dent et al., 1997;
   [221]Vassilatis et al., 1997; [222]Horoszok et al., 2001; [223]Hibbs
   and Gouaux, 2011), H. contortus ([224]Delany et al., 1998;
   [225]Forrester et al., 1999; [226]Cheeseman et al., 2001; [227]Portillo
   et al., 2003; [228]Yates et al., 2003; [229]McCavera et al., 2009) and
   Cooperia oncophora ([230]Njue and Prichard, 2004). However, GluCls have
   not been well characterized in T. canis. The T. canis genome encodes
   seven GluCls and these were expressed in the larval stages in this
   study. We clarified the nomenclature of GluCls using phylogenetic
   analysis ([231]Fig. 6A–[232]Table 2) as recommended by ([233]Beech et
   al., 2010). While GluCls were not among the DEGs in our initial
   analysis, upon further exploration we found that Tca-glc-2 was
   downregulated in ivermectin-treated larvae compared to the
   moxidectin-treated group ([234]Fig. 8). Additionally, Tca-glc-3 and
   Tca-avr-14 were downregulated in larvae exposed to ivermectin compared
   to controls ([235]Fig. 7, [236]Fig. 8). This finding is consistent with
   earlier work in Cooperia oncophora and Ostertagia ostertagi, where
   resistant isolates showed reduced avr-14B transcription compared to
   susceptible populations ([237]El-Abdellati et al., 2011). In H.
   contortus, polymorphisms and altered expression of avr-14 have also
   been implicated in resistance, though with considerable variability
   across isolates ([238]Glendinning et al., 2011). Experimental work in
   C. elegans has demonstrated that resistance requires simultaneous
   disruption of multiple GluCl α subunits, including avr-14, avr-15, and
   glc-1 ([239]Dent et al., 2000). Moreover, similar transcriptomic
   analysis in ivermectin-resistant strains of H. contortus from China
   showed that Hco-glc-5 was downregulated compared to susceptible strains
   from Australia ([240]Tuersong et al., 2022). However, studies have
   shown that Hco-glc-5 has no evidence of linkage to an ivermectin
   resistance phenotype in United Kingdom H. contortus isolates
   ([241]Laing et al., 2016; [242]Rezansoff et al., 2016). Further
   research is needed to clarify the specific roles of GluCls in ML
   resistance.

   ML resistance and tolerance in parasitic nematodes have been associated
   with P-gps (Synonym: ABCB1) ([243]Xu et al., 1998; [244]Kerboeuf et
   al., 2003; [245]Kerboeuf and Guegnard, 2011; [246]Peachey et al.,
   2017). P-gps have two ATP binding domains and two transmembrane domains
   ([247]Chen et al., 1986). Since P-gps are monocistronic, each fully
   annotated gene corresponds to a single protein. Other ABCB transporters
   include subclasses A, C, D, E, F, G, and H. The T. canis genome encodes
   thirteen P-gps and six other ABCB genes. Of these six additional ABCB
   genes, five were suitable for phylogenetic analysis after excluding
   Tca-abcb1 due to its short sequence. We clarified the nomenclature of
   these five genes using phylogenetic analysis ([248]Fig. 6B–[249]Table
   3). These genes were identified as half ABC transporters (haf).
   Tca-haf-10 did not cluster with any known C. elegans haf genes, and
   since no corresponding haf-10 gene exists in C. elegans, we designated
   it as Tca-haf-10. We retained the name Tca-abcb7, as this nomenclature
   is consistent with the orthologous gene in C. elegans. We did not
   identify any orthologs to Cel-haf-1, Cel-haf- 2, Cel-haf-3, Cel-haf-7,
   or Cel-haf-8 in T. canis. We assessed the expression level of the
   eighteen genes, including the thirteen named P-gp genes ([250]Jesudoss
   Chelladurai et al., 2023) and the five ABCB genes ([251]Table 3). While
   P-gps (ABCB1s) were not among the DEGs, Tca-haf-9 was upregulated in
   the moxidectin treatment group ([252]Supplemental Table 4A). Upon
   further analysis, differences in the expression of eight ABCB
   transporter genes were observed. In response to MLs, Tca-Pgp-11.2,
   Tca-Pgp-2, Tca-abcb7, Tca-haf-10, and Tca-Pgp-13.2, were downregulated,
   while Tca-Pgp-11.3, Tca-Pgp-16.3, and Tca-haf-9 were upregulated
   ([253]Fig. 7, [254]Fig. 8). This is in agreement with previous findings
   that T. canis larvae exposed to ivermectin downregulate the expression
   of Tca-Pgp-13.1, while exposure to milbemycin upregulates Tca-Pgp-16.2
   and Tca-Pgp-16.3. ([255]Jesudoss Chelladurai et al., 2023).

   Beyond GluCl and P-gp genes of interest, additional pathways impacted
   by MLs include those involved in the regulation of transcription,
   energy production, structural and functional roles, essential
   physiological functions like reproduction, expression of genes
   responsible for excretory/secretory (ES) molecules, and other processes
   of potential biological significance.

   This study suggests that MLs influence gene transcription by modulating
   mechanisms that control transcriptional regulation. In this study,
   ivermectin was observed to upregulate a mediator of RNA polymerase II
   transcription subunit 13, Tca-let-19 (Tcan_05892). This gene is part of
   the mediator complex that governs gene expression by influencing
   transcriptional initiation and regulation ([256]Poss et al., 2013).
   Similar findings in osteosarcoma models suggest that this mediator
   complex plays a critical role in promoting tumor growth and cellular
   stress responses ([257]Yu et al., 2014). In addition, ivermectin
   upregulated the putative U5 small nuclear ribonucleoprotein helicase,
   Tca-Y46G5.4 (Tcan_05430), which is involved in pre-mRNA splicing and
   post-transcriptional RNA modification. Moxidectin, on the other hand,
   downregulated the putative zinc finger protein F56D1.1 (Tcan_16300), a
   protein known to bind DNA, RNA, and proteins, playing a vital role in
   gene regulation, DNA repair, and other essential cellular processes. GO
   enrichment analysis provided further evidence that several
   transcription-related processes were significantly affected, including
   the upregulation of pathways associated with nucleosome positioning,
   RNA destabilization, ribosomal activity, and translational initiation
   in larvae exposed to MLs ([258]Fig. 3, [259]Fig. 4B). These demonstrate
   the impact MLs have either directly or indirectly on transcriptional
   regulation.

   Our findings indicate ML exposure may impact energy metabolism by
   modulating the expression of key metabolic enzymes. Ivermectin and
   moxidectin upregulated mitochondrial glutamate dehydrogenase 2
   (Tca-GLUD2, Tcan_03811), an enzyme studied in the parasitic nematodes
   Teladorsagia circumcincta and H. contortus that is crucial for nitrogen
   and energy metabolism, converting glutamate to α-ketoglutarate and
   ammonia ([260]Skuce et al., 1999; [261]Umair et al., 2011). In
   contrast, ivermectin downregulated methylmalonyl-CoA epimerase,
   mitochondrial (Tca-Mcee, Tcan_12863), an enzyme involved in the
   metabolism of fatty acids and amino acids in C. elegans ([262]Kuhnl et
   al., 2005). Ivermectin also downregulated
   amidophosphoribosyltransferase (Tca-Prat, Tcan_18524), which catalyzes
   the first step in purine biosynthesis and has been shown to be
   differentially expressed in Wolbachia during development of Brugia
   malayi male and female adults ([263]Grote et al., 2017). Moxidectin
   downregulated other enzymes involved in amino acid and vitamin
   metabolism, including a lysine histidine transporter 1 (Tca-LHT1,
   Tcan_17795), a membrane protein responsible for the transport of lysine
   and histidine; a mitochondrial alpha-aminoadipic semialdehyde synthase,
   (Tca-AASS, Tcan_18345), shown to be involved in lysine degradation in
   C. elegans ([264]Zhou et al., 2019); and a mitochondrial methylmalonic
   aciduria and homocystinuria type D-like protein, (Tca-Mmadhc,
   Tcan_06175), which plays a role in Vitamin B12 metabolism. GO
   enrichment analysis further revealed that DEGs in response to MLs were
   related to metabolic pathways, including upregulation of those
   associated with NADH dehydrogenase activity and the respiratory chain
   ([265]Fig. 3, [266]Fig. 4B), indicating MLs may affect T. canis’
   ability to metabolize energy efficiently.

   MLs appear to play a role in body structure and function by modulating
   the transcription of key genes. Ivermectin upregulated nuclear
   anchorage protein 1 (Tca-anc-1, Tcan_00659), which is essential for
   maintaining nucleus-cytoskeleton interactions and stabilizing neuron
   positions in C. elegans ([267]Johnson and Kramer, 2012). Moxidectin
   similarly upregulated protocadherin-16 (Tca-F32A5.4, Tcan_07881), a
   protein likely involved in cell-cell adhesion and neuronal development,
   contributing to tissue morphogenesis and homeostasis. In contrast,
   moxidectin downregulated cuticlin-1 (Tca-cut-1, Tcan_11029), described
   in the parasitic nematode Ascaris suum as a key structural protein of
   the epicuticle that is rich in proline and alanine ([268]Fujimoto and
   Kanaya, 1973; [269]Bisoffi et al., 1996). In the parasitic nematode B.
   malayi, downregulation of cuticle collagen genes in response to
   ivermectin has been studied ([270]Ballesteros et al., 2016).
   Furthermore, in ivermectin-resistant H. contortus, studies have shown
   that cuticle collagen genes are downregulated when compared to
   susceptible isolates ([271]Tuersong et al., 2022). Additionally,
   moxidectin downregulated TRIO and F-actin-binding protein (Tca-Triobp,
   Tcan_10839). In C. elegans, UNC-73, a TRIO homolog, is involved in
   regulating cytoskeletal and neuronal development ([272]Steven et al.,
   2005). Finally, the differential expression of genes corresponding to
   the GO enrichment pathways for the troponin complex, sarcoplasmic
   reticulum, and myofilament formation ([273]Fig. 3B) suggest MLs may
   disrupt normal muscular processes.

   MLs may impact the transcription of genes involved in physiological
   functions, including reproduction. Ivermectin downregulated chondroitin
   proteoglycan 4 (Tca-cpg-4, Tcan_10605), which may play a role in
   reproductive processes, similar to chondroitin proteoglycan 2
   (Tca-cpg-2); Tca-cpg-2 in adult female T. canis is crucial for
   oogenesis and embryogenesis and is enriched in germline tissues
   ([274]Ma et al., 2018). Additionally, ivermectin downregulated
   nanos-like protein 2 (Tca-Nanos2, Tcan_18885), which is involved in
   germ cell development. GO enrichment analysis showed that MLs
   downregulated pathways involved in cellular responses to hormone
   stimuli and longevity-regulating processes which could suggest that MLs
   interfere with essential biological functions ([275]Fig. 4B).

   Our results suggest that MLs impact the transcription of genes that
   encode ES and other molecules crucial for host-parasite interactions.
   Ivermectin downregulated a metalloproteinase inhibitor 2 (Tca-TIMP2,
   Tcan_03888). In Ancylostoma caninum, the metalloproteinase inhibitor
   AcTMP is speculated to regulate host matrix metalloproteinases or host
   neutrophil collagenase activity to influence tissue repair and
   inflammation at the intestinal mucosa attachment site ([276]Knox,
   2007). Moxidectin, on the other hand, upregulated multiple protective
   factors, including a trypsin inhibitor (Tcan_08038) and an aspartyl
   protease inhibitor (Tca-API-3, Tcan_02586), which may protect against
   host-specific proteolytic enzymes and play a role in the host immune
   response and infection, as seen in A. suum ([277]Hawley and Peanasky,
   1992), H. contortus ([278]Li et al., 2017), and Ostertagia ostertagi
   ([279]De Maere et al., 2005). Additionally, moxidectin upregulates a
   group XV phospholipase A2 (Tca-PLA2G15, Tcan_15544). In parasitic
   nematodes, phospholipase A2 (PLA2) plays crucial roles in pathogenesis
   by contributing to host tissue penetration, modulating host immune
   responses, and facilitating nutrient acquisition ([280]Belaunzarán,
   2023). These enzymes disrupt host cell membranes, release signaling
   molecules that manipulate host immune pathways, and help the parasites
   maintain cellular integrity and repair, thereby enhancing their
   survival and ability to infect the host ([281]Belaunzarán, 2023).
   Finally, moxidectin downregulated arginase (Tca-rocF, Tcan_06332),
   which normally competes with nitric oxide synthase for L-arginine,
   reducing nitric oxide production and aiding immune evasion [82]. The
   differential expression of genes responsible for the production of ES
   molecules may indicate that MLs alter the ability of T. canis to
   interact within the host.

   This study indicated that MLs influence other genes in T. canis.
   Ivermectin downregulated a putative oxidoreductase, Tca-dhs-27
   (Tcan_13635), an enzyme involved in redox reactions that is crucial for
   maintaining cellular redox homeostasis. Interestingly, a similar
   oxidoreductase is upregulated in Anisakis simplex larvae in response to
   ivermectin ([282]Polak et al., 2020).

   Ivermectin also downregulated UDP-glucuronosyltransferase 2C1,
   Tca-UGT2C1 (Tcan_05445), an enzyme for detoxifying both endogenous and
   exogenous compounds. UDP-glucuronosyltransferases facilitate the
   glucuronidation process, which makes drugs and xenobiotics more
   water-soluble, enhancing their excretion ([283]Dimunova et al., 2022).
   The involvement of these genes in T. canis larvae, and across life
   cycle stages, needs further investigation.

   MLs appear to affect the transcription of several unnamed genes, which
   have unknown functions (distinct from pseudogenes). Ivermectin
   upregulated Tcan_06163, identified as an ortholog of C. elegans
   T20G5.12, though its function remains unknown. Additionally, ivermectin
   upregulated Tcan_00335, a SET domain family protein. In C. elegans, SET
   domain proteins play a crucial role in regulating gene expression
   through histone methylation, influencing key processes such as
   development, reproduction, chromatin remodeling, and lifespan control
   ([284]Terranova et al., 2002; [285]Ni et al., 2012; [286]Engert et al.,
   2018). Ivermectin downregulated Tcan_14498, a domain-containing
   protein. In C. elegans, this domain has no known functions (pfam
   PF01682) but has been associated with immunoglobulin and fibronectin
   domains, and lipases ([287]Consortium, 2024). Moxidectin upregulated
   Tcan_00043, which was identified as either an ortholog of the V-type
   proton ATPase subunit A in D. immitis or Bm33 in B. malayi. V-type
   ATPases are involved in acidifying intracellular endosomes, lysosomes,
   and Golgi-derived vesicles and hyperpolarizing membranes through ATP
   hydrolysis to produce a proton gradient ([288]Abaza, 2022). Bm33 is an
   aspartyl protease inhibitor thought to play a role in immune evasion by
   inhibiting host proteases ([289]Maizels et al., 2001). Moxidectin also
   upregulated Tcan_00055, identified as an ortholog to a regulator of
   rDNA transcription protein 15 in Trichinella papuae ([290]Korhonen et
   al., 2016). Hypothetical proteins without identifiable orthologs were
   also observed. The lack of functional annotations for these genes
   limits our understanding of their precise roles within T. canis and
   therefore limits our ability to predict the impact of MLs.

   Some limitations of this study include the use of a small sample size
   (n = 3 replicates per treatment) and only a 12-h drug exposure period,
   which may not fully capture the range of biological responses to ML
   treatments. Furthermore, variability in the ivermectin-treated group is
   notable, with one replicate showing lower similarity to the others
   (intragroup correlation = 0.71) compared to higher correlations among
   the remaining samples. Additionally, the study was conducted in vitro
   using a single 10 μM dose which limits the generalizability of the
   findings to real-world scenarios where drug exposure may vary. While
   these results provide insights into ML responses in T. canis hatched
   larvae, extrapolation to hypobiotic larvae should be interpreted with
   caution. Another limitation is the absence of a well-annotated T. canis
   chromosome-level genome assembly, making accurate mapping more
   challenging. Improving the reference genome and its annotation would
   greatly enhance the resolution of transcriptomic analyses. Future
   research should investigate the effects of multiple doses over extended
   periods, explore alternative drugs, and incorporate multi-omics
   approaches, such as proteomics, to provide a more comprehensive
   understanding of the biological mechanisms involved. Additionally,
   mouse models could be employed to study and compare transcriptomic
   results from untreated and ML-treated larvae in vivo to in vitro
   studies, facilitating the interpretation of the impact of MLs on
   somatic L3s within a host.

5. Concluding statement

   In conclusion, this study provides valuable insights into how MLs like
   ivermectin and moxidectin impact the transcription of various genes
   associated with essential biological processes in T. canis hatched L3s.
   In L3s treated in vitro with ivermectin or moxidectin, we identified
   608 and 1413 DEGs, respectively, of the 18,596 genes in the genome.
   Among the transcriptome were genes known to be associated with ML
   resistance, including GluCls and P-gps, which remain of significant
   interest as potential targets for therapeutic interventions. We also
   identified DEGs of potential interest that could confer ML tolerance in
   T. canis L3s, which could serve as valuable candidates for future
   research and exploration. Future studies should employ the use of
   multi-omics approaches to capture the full spectrum of molecular
   changes induced by MLs. Additionally, investigating the effects of MLs
   at multiple doses over extended periods will help determine long-term
   impacts on T. canis larvae. Researchers should also consider exploring
   the impacts of alternative anthelmintic compounds to identify
   treatments that may overcome tolerance mechanisms associated with MLs.
   Lastly, comparing in vivo larvae to in vitro larvae studies will help
   clarify the specific gene expression changes and physiological effects
   induced by MLs within a host environment. Together, these approaches
   will provide a comprehensive understanding of MLs’ effects on T. canis.

CRediT authorship contribution statement

   Theresa A. Quintana: Writing – review & editing, Writing – original
   draft, Visualization, Validation, Software, Methodology, Investigation,
   Formal analysis, Data curation, Conceptualization. Matthew T. Brewer:
   Writing – review & editing, Conceptualization. Jeba R.J. Jesudoss
   Chelladurai: Writing – review & editing, Writing – original draft,
   Visualization, Validation, Supervision, Software, Project
   administration, Methodology, Investigation, Funding acquisition, Formal
   analysis, Data curation, Conceptualization.

Ethics statement

   All protocols were approved by the Institutional Biosafety Committee at
   Kansas State University (IBC-1516-VCS)

Funding

   This project was funded by the NIH Grant CEZID 5P20GM130448 and
   start-up funds provided to JJC by the College of Veterinary Medicine
   Kansas State University.

Conflicts of interest statement

   I, Jeba R J Jesudoss Chelladurai, on behalf of the authors of this
   manuscript, hereby declare that all the authors who have contributed to
   this manuscript have no conflicts of interest to disclose in relation
   to this submission. I confirm that there are no financial,
   professional, or personal relationships that could be perceived as
   influencing the content or outcomes of the research presented in the
   manuscript titled "Transcriptional responses to in vitro macrocyclic
   lactone exposure in Toxocara canis larvae using RNA-seq."

   I also affirm that all contributions to this work, including research,
   data analysis, and writing, were carried out in an unbiased and
   transparent manner.

Acknowledgements