P-glycoproteins and amphids: a two-stage trajectory of ivermectin resistance in Caenorhabditis elegans

TL;DR

This study shows how C. elegans develops resistance to ivermectin in two stages, involving early drug efflux and later changes in sensory structures.

Contribution

The paper identifies a two-stage mechanism of ivermectin resistance involving P-glycoproteins and amphid defects in C. elegans.

Findings

Early ivermectin exposure triggers transient P-glycoprotein overexpression.

Long-term resistance is linked to ciliary/neuronal pathway remodeling and a Dyf phenotype.

Three independent lineages converged on similar resistance mechanisms despite genetic differences.

Abstract

Ivermectin (IVM) is a cornerstone of nematode control. However, its efficacy is increasingly compromised by emerging resistance in target parasites. P-glycoproteins (Pgps), which mediate drug efflux, and amphid neuron defects, characterized by a dye-filling defective (Dyf) phenotype, may both contribute to IVM resistance in Caenorhabditis elegans, but their respective roles and potential connection remain to be clarified. Here, we investigated these mechanisms in three independently evolved Caenorhabditis elegans lineages exposed to stepwise IVM selection. These populations, although distinct, converged over time to a comparable and stabilized high level of IVM resistance. Endpoint transcriptomes revealed fluctuations in the expression of several pgp, overexpressed in one lineage only, suggesting that sustained pgp upregulation is not a conserved feature of stable resistance. Consistently, neither short-term exposure to IVM nor long-term drug pressure selection elicited strong transcriptional induction. Nevertheless, we uncovered a transient surge of several pgp transcripts during the early selection process of the worms exposed to low doses of IVM. We speculate this relates to an early pgp-dependent tolerance mechanism in the dynamic adaptation program to IVM. Beyond an initial pgp-based response, our expanded RNA-seq analysis across lineages revealed a conserved enrichment for ciliary/neuronal pathways including genes expressed in amphid/phasmid, indicating remodeling of the chemosensory/ciliary network as resistance intensifies. All resistant lines also converged on the Dyf phenotype that emerged during selection at stages coinciding with sharp gains in IVM tolerance, implicating amphid structure/function in the establishment of durable resistance. Together, these data reconcile heterogeneous literature on PGPs by decoupling transient, context-dependent pgp activation from long-term maintenance of resistance, and are consistent with previous work showing ciliary pathways as common denominators of IVM resistance. We propose a two-phase associative model in which early PGP-mediated tolerance precedes the emergence of ciliary/amphidial signatures and a Dyf phenotype during selection that durably reduces susceptibility.

Image 1

•Comparison of three independently evolved ivermectin-resistant C. elegans lines.•Two-stage trajectory of ivermectin resistance acquisition.•Early, transient pgp overexpression at low ivermectin doses.•Cross-lineage RNA-seq highlights conserved ciliary/neuronal signatures.•Dyf phenotype co-occurs with sharp increases in ivermectin tolerance.

Comparison of three independently evolved ivermectin-resistant C. elegans lines.

Two-stage trajectory of ivermectin resistance acquisition.

Early, transient pgp overexpression at low ivermectin doses.

Cross-lineage RNA-seq highlights conserved ciliary/neuronal signatures.

Dyf phenotype co-occurs with sharp increases in ivermectin tolerance.

Figures

Fig. 4Temporal transcription profile of *pgp* expression following withdrawal of IVM pressure.(A) The IVR10-2009 strain was maintained for up to 6 months either under IVM pressure (10 ng/mL IVM; denoted IVR10(+)) or without selection (0 ng/mL IVM; denoted IVR10(−)). Following this initial phase, defined as time point M0, each condition was either continued or switched to the opposite treatment (±IVM) for an additional 6 months, resulting in four selection pressure histories: IVR10(+,+); IVR10(+,-); IVR10(-,+) and IVR10(-,-). *pgp* mRNA levels were measured by real-time quantitative PCR (RT-qPC

Amphid dye uptake.Adult worms from the selection process of IVR10-2014 (A) or IVR10-2022 (B) and corresponding to time points: early stage of the 2 ng/mL plateau (2a), end of the 2 ng/mL plateau (2b), 6 ng/mL (6), and the end of the last plateau at 10 ng/ml (10) for IVR10-2014 and 11 ng/mL (11) for IVR10-2022 ([Fig. 5](#fig5)A and B) were incubated with the red fluorescent lipophilic dye Dil.Fig. 7

Full text

Introduction

1

Several billion people live in areas at risk of infection by parasitic nematodes that cause chronic and debilitating diseases. Nematodes also significantly impact the health and welfare of companion animals, livestock and cause significant economic losses in agriculture. Their control relies on the use of anthelmintics, mainly from the macrocyclic lactone (ML) family, of which ivermectin (IVM) is the most widely used (Martin et al., 2021). MLs act by targeting glutamate-gated chloride channels in nematodes, leading to neuronal hyperpolarization and paralysis (Dent et al., 2000). Reduced efficacy and resistance signals have been reported across a broad range of nematode infections in human and veterinary medicine, including gastrointestinal strongyles of ruminants (e.g., Haemonchus contortus), filarial nematodes affecting humans (e.g., Onchocerca volvulus), and canine heartworm (Dirofilaria immitis), among others. Accordingly**,** increasing treatment failures across these species raise serious concerns about the sustainability of current control strategies (Kotze and Prichard, 2016; Martin et al., 2021; Prichard, 2021; Specht and Keiser, 2023). Therefore, it is urgent to understand how ML resistance emerges and stabilizes across nematode lineages.

The nematode Caenorhabditis elegans is a widely used nematode model for investigating anthelmintic resistance due to its genetic tractability, short life cycle, and high susceptibility to a broad range of antiparasitic compounds, including IVM. Its well-annotated genome, self-fertilizing reproduction, and small effective population size make it particularly suitable for controlled studies of transcriptional reprogramming and stochastic evolutionary dynamics under defined drug selection pressures. Moreover, functional conservation with parasitic nematodes allows C. elegans to serve as a powerful proxy for identifying anthelmintic resistance-associated genes (Hahnel et al., 2020).

Mutations in three glutamate-gated chloride channel (GluCl) genes (glc-1, avr-14, and avr-15) is required for high levels of resistance to IVM in C. elegans (Dent et al., 2000). While this finding validate the pharmacological target and underscore the utility of C. elegans for mechanism dissection, their field relevance appears limited. Indeed, mutations in GluCls have not been consistently associated with H. contortus resistant populations (Doyle et al., 2022; Kotze et al., 2014; Williamson et al., 2011), supporting the view that other non-targeted mechanisms predominate in natural parasite populations.

The molecular basis of ML resistance in nematodes is currently recognized as multigenic. Across species, it is thought that ML resistance involves coordinated changes in drug disposition and sensory/structural pathways, including increased xenobiotic efflux via ATP-binding cassette (ABC) transporters (notably P-glycoproteins, PGPs), altered biotransformation (e.g., cytochrome P450–mediated phase I metabolism and phase II conjugation), and ciliated sensory neurons defects that influence exogenous molecules uptake (Freeman et al., 2003; Laing et al., 2022; Lespine et al., 2024; Urdaneta-Marquez et al., 2014). Among non-target mechanisms, PGPs have received considerable attention. In mammals, these efflux pumps actively export a broad range of endogenous and exogenous compounds including IVM, out of cells, thereby altering homeostasis, reducing toxicity and limiting drug efficacy. PGP structures are highly conserved across kingdom supporting their drug efflux function in nematodes (Jin et al., 2012). In C. elegans and parasitic nematodes, several mutation in pgp genes have been associated to IVM resistance, and pgp have been shown to be overexpressed in resistant populations (Lespine et al., 2024). Across parasitic nematodes, pgp-9 is implicated in avermectin interaction and resistance (e.g., Haemonchus contortus; Teladorsagia circumcincta), while pgp-11 carries resistance-linked polymorphisms and functionally reduces IVM susceptibility when expressed transgenically (e.g., Dirofilaria immitis; Parascaris) (Curry et al., 2022; Dicker et al., 2011; Godoy et al., 2016; Janssen et al., 2015; Mani et al., 2016; Turnbull et al., 2018). In C. elegans, regulation of pgp-1, pgp-3, pgp-6, pgp-9, and *pgp-*13 by the conserved nuclear hormone receptor NHR-8 has been proposed to modulate IVM susceptibility. Accordingly, loss of nhr-8 sensitizes worms to IVM (Ménez et al., 2019). These results suggest that IVM exposure induces pgp transcriptional overexpression via gene transcriptional regulation, contributing to IVM tolerance. However, findings on pgp expression in resistant nematodes vary qualitatively and quantitatively and the dynamics of pgp regulation in resistance remain poorly defined (Lespine et al., 2024).

In parallel to efflux-based mechanisms, a growing body of literature shows that disruption of amphid cilia via mutations in intraflagellar transport genes is tightly linked to increased resistance to ML in C. elegans (Brinzer et al., 2024; Page, 2018). 25 years ago, Dent et al. identified dye-filling defective (Dyf) mutants in a genetic screen for resistance to 10 ng/ml of IVM in an avr-15 mutant background (Dent et al., 2000). The Dyf phenotype is defined by the inability of worms to take up environmental molecules, such as the lipophilic fluorescent dye Dil. It results from structural defects in ciliated sensory neurons, particularly the amphids. Recent studies have confirmed a strong association between the Dyf phenotype and IVM resistance (Brinzer et al., 2024; Page, 2018).

Despite these findings, some questions remain such as the magnitude and context-dependence of PGPs’ contribution, and how detoxification networks integrate with defects in ciliated sensory neurons. Clarifying these unknowns requires tractable models enabling controlled selection and time-resolved molecular phenotyping.

To answer these questions, several ML-resistant lines of C. elegans were generated independently by selection upon stepwise exposure to IVM (James and Davey, 2009). Transcriptional analyses and Dil staining were performed on these strains to compare pgp gene expression levels and amphids function along the selection process.

Our results show that pgp upregulation is lineage-specific and transient during resistance emergence, but is not always maintained once resistance has stabilized. Furthermore, pgp are not strongly induced by acute or sustained IVM exposure. By expanding our analysis of resistance-associated genes through RNA sequencing, we found an enrichment of deregulated genes in ciliary/neuronal pathways common to all resistant strains. Accordingly, all lines exhibit a Dyf phenotype that emerged when resistance to high doses of IVM developed, providing a common signature in the evolution of IVM resistance.

Materials and methods

2

C. elegans nematode strains and cultivation conditions

2.1

Wild-type C. elegans Bristol strain N2 was obtained from the Caenorhabditis elegans Genetics Center (CGC; University of Minnesota, Minneapolis, MN, USA). The IVR10-2009 and IVR-10-2014 strains, selected from the wild-type strain with IVM and phenotypically resistant to IVM, were respectively kindly provided by C. E. James (James and Davey, 2009) and previously selected in our laboratory (Ménez et al., 2016). C. elegans strains were cultured and handled according to the procedures previously described (unless otherwise stated) (Brenner, 1974). Briefly, nematodes were cultured at 21° on Nematode Growth Medium (NGM) agar plates (1.7% bacto agar, 0.2% bactopeptone, 50 mM NaCl, 5 mg/L Cholesterol, 1 mM CaCl2, 1 mM MgSO4, and 25 mM KPO4 Buffer) seeded with Escherichia coli OP50 strain as a food source. IVM-containing NGM plates were prepared as follows: stock solutions of IVM in DMSO were diluted in NGM to the appropriate concentration and at a final DMSO concentration of 0.1% before pouring plates. The IVR10-strains were cultured on NGM plates containing 11.4 nM (10 ng/ml) of IVM, unless otherwise stated. Nematodes were synchronized through egg preparation with sodium hypochlorite. Briefly, 5 to 15 gravid adults were picked on a NGM plate without IVM. After one generation, asynchronous populations, consisting mainly of gravid adults and eggs were collected by washing the surface of the NGM plates with M9 buffer (3 g KH2PO4, 6 g Na2HPO4, 5 g NaCl, and 0.25 g MgSO4 in 1 L of water) and centrifuged at 1200×g for 30 s. All larval stages except eggs were lysed with a bleaching mixture (5 M NaOH and 1% sodium hypochlorite) during 5 min. After three washes with M9 buffer, C. elegans eggs were hatched overnight at 21 °C in M9 solution without bacteria to obtain a synchronized first-stage larval (L1) population.

Stepwise exposure of C. elegans N2 Bristol to IVM

2.2

Culture conditions for the development of ML-resistant C. elegans strains following stepwise exposure to IVM were previously described (James and Davey, 2009; Ménez et al., 2016). Briefly, at week 0, hundreds of wild-type worms were transferred to NGM plates containing 0.5 ng/ml IVM for IVR-2022. Twice weekly, worms were transferred to fresh NGM plates. Once the worms survived and reproduced, they were transferred to plates containing higher doses of IVM. The IVM concentrations used to generate IVM-selected strains from wild-type N2 Bristol were 0.5, 1, 2, 3, 4, 5, 6, 7, 8, 9 and 10 ng/ml. After approximately 27 to 44 weeks, wild-type worms were able to survive and to be maintained on 10 ng/ml of IVM.

Larval development assay (LDA) on C. elegans strains

2.3

LDA measures the potency of anthelmintics to inhibit the development of C. elegans from L1 to the young adult stage. Approximately 30 synchronized L1-stage larvae were added to each well of a 12-well plate containing NGM medium supplemented with increasing concentrations of compound of interest and seeded with OP50 bacteria. DMSO vehicle was used as control at a maximal concentration of 0.3%, a level at which no harmful effect on C. elegans were observed. Plates were then incubated at 21° in the dark for 52–55 h, during which the L1 larvae in the negative control were developed into late L4/young adult worms. L1, L2 and L3 stages were considered inhibited in their development, while late L4 and young adult worms were classified as developed. Development was calculated as the percentage of worms reaching the late L4 or young adult stage in the presence of IVM, normalized to the untreated control. Each concentration was set-up in triplicate, and all experiments were independently replicated at least three times. A sigmoidal dose-response curve with variable slope was fitted using Prism 8 (GraphPad, SanDiego, USA), allowing calculation of the effective concentration for a 50% effect (EC_50_). Statistical significance was assessed using Student's t-test.

Transcriptome analysis

2.4

RNA preparation

2.4.1

Approximatively 2000 synchronized L1 larvae were added in NGM plates containing or not variable concentrations of IVM and seeded with OP50 bacteria. After 40-45 h, L4-stage worms were collected using a bovine serum albumin coated 1000 μl pipette tip. The worms were washed three times with M9 buffer. Worms were then resuspended in 600 μl of RLT buffer, supplied in the RNeasy Mini Kit (Qiagen, Hilden, Germany), with 12 μl of 2 M DTT and flash-frozen in liquid nitrogen. The worms were homogenized using a FastPrep-24 (MP Biomedicals, Irvine, USA) with lysis matrix D; three runs of 10 s at 6 m/s were performed, with 2 min of cooling on ice between runs. Subsequently, RNA was extracted using RNeasy Mini Kit (Qiagen, Hilden, Germany). RNA quality for RNA-seq was assessed using the TapeStation 4200 (Agilent, Santa Clara, USA). For RT-qPCR, total RNA was quantified with a NanoDrop One spectrophotometer (NanoDrop Technologies Inc., Wilmington, USA), and RNA purity was confirmed by measuring the A260 nm/A280 nm ratio, which consistently ranged between 1.8 and 2.0.

Gene expression studies by RNA-seq

2.4.2

For each of the samples, RNA-seq libraries were constructed from 1 μg of total RNA at the GeT-TRiX facility (GénoToul, Génopole Toulouse Midi-Pyrénées) using Illumina® Stranded mRNA Prep kit (Illumina, San Diego, USA) following the manufacturer's instructions; with an adaption to produce library sizes compatibles with paired-end 150 bp read length sequencing. Briefly, mRNAs were captured, fragmented, and converted to double-stranded cDNA, followed by ligation of sequencing adaptors. Libraries were amplified with eleven cycles of PCR.

The yield and quality of the libraries were assessed using D1000 ScreenTapes on TapeStation 4200 (Agilent Technologies, Santa Clara, USA). The libraries were subsequently pooled at equimolar concentrations and transferred to GeT-PlaGe facility (GénoToul, Génopole Toulouse Midi-Pyrénées) for sequencing. The pooled library was loaded onto a single sequencing lane on an Illumina NovaSeq 6000 S4 v1.5 flowcell and sequenced using a 2 × 150 bp paired-end sequencing mode.

Bioinformatics processing of the sequencing data was performed by the GeT-TRiX facility (GénoToul, Génopole Toulouse Midi-Pyrénées) using the computing resources of Genotoul Bioinformatics platform (Bioinfo Genotoul, Toulouse Midi-Pyrenees). The analysis pipeline was executed with Nextflow v22.12.0-edge (Di Tommaso et al., 2017) and processed using nf-core/rnaseq v3.12.0 (Patel et al., 2023) from the nf-core collection of workflows (Ewels et al., 2020). Reads were aligned to C. elegans reference genome PRJNA13758WBPS18.

Biostatistics analyses were performed using R v4.2.2 (R Core Team, 2022), the edgeR package (Robinson et al., 2010), and Bioconductor packages (Huber et al., 2015). Briefly, the raw count table was filtered to remove undetected genes using filterByExpr function with min.count = 10. Normalization factors were calculated using TMM method. Common, trended, and tagwise negative binomial dispersions were estimated with the robust mode of estimateDisp function. A negative binomial generalized log-linear model was fitted to the count data with the quasi-likelihood (QL) methods, with glmQLFit function to moderate the genewise QL dispersion (McCarthy et al., 2012). Pair-wise comparisons between biological conditions were performed using specific contrasts with glmQLFtest function to identify differentially expressed genes. Multiple testing correction was applied using the Benjamini-Hochberg procedure (Benjamini and Hochberg, 1995) to control the false discovery rate (FDR). Genes with an FDR ≤0.05 were considered to be differentially expressed between conditions. RNAseq data are publicly available in the Gene Expression Omnibus (GEO) database under accession number [GSE314492](GSE314492).

RNA-seq-based functional enrichment analysis

2.4.3

Following differential transcriptional expression analysis (as described in the corresponding section), significantly upregulated and downregulated genes (FDR ≤0.05) were identified separately for each IVM-resistant C. elegans strain. Functional enrichment analysis was performed using ShinyGO v0.82, with gene identifiers mapped to Gene Ontology (GO) terms and KEGG pathways using annotation data from WormBase and the STRING database.

Quantification of mRNA transcripts by RT-qPCR

2.4.4

cDNA was synthesized from 2 μg of RNA using the Maxima H Minus cDNA Synthesis Kit (Thermo Fisher, Waltham, MA, USA), following the manufacturer's instructions. RT-qPCR was performed on a QuantStudio 5 Detection System (Applied Biosystems, Courtaboeuf, France) with Power SYBR Green Master Mix 2X (Applied Biosystems).

Gene-specific primers (see Ménez et al., 2019) and primers for tba-1 (forward: ATCGATTTTTGTAGATCTTGAGCCA; reverse: TCCAGTGCGGATCTCATCAAC) were designed using Primer Express v2.0 software (Applied Biosystems), based on the C. elegans genome sequence available at www.wormbase.org. Primer sequences were synthesized by Invitrogen (Cergy-Pontoise, France), and their specificity was confirmed using the NCBI BLAST tool. Raw fluorescence data were analyzed with LinRegPCR, which performs baseline correction and determines amplification efficiency from the exponential phase. Individual reaction efficiencies were incorporated into relative quantification. Mean starting quantities (N_0_ values) were normalized to the reference gene cdc-42 or tba-1. Statistical analysis was conducted with GraphPad Prism (San Diego, USA) using Student's t-test.

IVM exposure

2.5

Short-term IVM exposure and RT-qPCR analysis

2.5.1

Approximately 2000 synchronized L1 larvae of N2 Bristol and IVR10-2009 strains were plated on NGM plates seeded with E. coli OP50 and incubated at 21 °C until L4 stage (40–45 h). Worms were washed with M9 buffer and transferred to NGM plates containing DMSO (control) or increasing concentrations of IVM, selected based on IVM EC_50_ for each strain and phenotypic thresholds (immobility, clustering). Concentrations ranged from 2.5 to 30 nM for N2 and 7.5 to 100 nM for IVR10-2009. After 6 h, worms were harvested for RNA extraction.

pgp transcript levels were measured by RT-qPCR, normalized to tba-1, and expressed as fold change relative to DMSO control (set to 1). Results are presented as mean ± S.D. from 3 to 4 independent experiments. Statistical significance was assessed using one-way ANOVA with Dunnett's test.

Long-term IVM selection and withdrawal experiments

2.5.2

To assess the effect of prolonged IVM selection or withdrawal on pgp transcriptional expression, IVR10-2009 was maintained for 6 months on NGM plates with (IVR10-2009+) or without (IVR10-2009−) 10 ng/mL IVM. Worms were then either kept under the same condition or switched, generating four IVR10 histories: +,+ (constant), +,− (withdrawn), −,+ (added), −,− (permanently removed).

Samples were collected at: 0, 3, 6 months, and weekly (W0–W4) or monthly (M0–M6). RNA was extracted from mixed-stage worms, and pgp expression was measured by RT-qPCR, normalized to cdc-42, and expressed relative to the baseline T0 (set to 1). Values represent one (Fig. S3) or three (Fig. 4) independent biological replicates. Two-way ANOVA with Dunnett's test was used.

C. elegans dye-filling assay (Dil staining of amphids)

2.6

To visualize the amphids, worms were synchronized at the adult stage. Worms were then incubated in a dye solution containing 10 ng/ml of DilC12(3) in M9 broth with gentle shaking for 2 h at 21 °C. After a recovery period of 2 h on NGM plates, worms were paralyzed by using levamisole (40 mM), and observed by using a Nikon SMZ800N microscope with a fluorescent filter for RFP, equipped with a VisiCam 5 Plus camera and analyzed by using waveimage software.

Results

3

Adaptation of C. elegans to IVM pressure through generations

3.1

To study the mechanisms associated with IVM resistance we used three independent C. elegans strains selected under increasing IVM pressure. As a reference resistant strain, IVR10-2009 was originally generated and characterized by James and Davey in 2009 in Australia (James and Davey, 2009). Following the same gradual exposure approach, we independently selected two other resistant strains in our laboratory. The IVR10-2014 line described in a previous study (Ménez et al., 2016), and IVR10-2022, a recently developed strain. We plotted the time-course of selection of worms exposed to increasing IVM concentration which reflected their progressive adaptation to the drug (Fig. 1A). Interestingly, all three strains showed progressive adaptation to higher IVM doses following comparable dynamics with successive plateaus but with subtle differences. A first plateau at 2 ng/ml was common to all lines followed by additional plateaus at 6–7 ng/ml observed only for IVR10-2009 and IVR10-2014 lines. These plateaus reflect the inability of the worm population to withstand higher IVM doses. After 39 to 48 weeks of drug pressure, the selection process was ended at a final 10-11 ng/ml plateau for all IVM-selected strains. At that stage, the worms were maintained at a constant IVM concentration of 10 ng/ml. All strains were able to survive and to reproduce in these conditions.Fig. 1Generation of IVM-resistant C. elegans strains and quantification of relative mRNA levels of pgp by RNA-seq.(A) Time-course of stepwise selection of IVM resistance from wild-type N2 Bristol C. elegans strain. Worms, cultured on Nematode Growth Medium (NGM) plates containing IVM, were transferred to NGM plates containing higher doses of IVM once they were able to survive and reproduce. The selection curves of the IVR10-2014 and IVR10-2022 lines reflect the acquisition of resistance in our laboratory, while the selection curve of the IVR10-2009 strain, kindly provided by C. E. James, displays the timeline described by James and Davey in 2009 (James and Davey, 2009). (B) Sensitivities of N2 Bristol and IVM-selected C. elegans strains to IVM assessed by larval development assay. Values represent the percentage of L1 larvae reaching the young adult stage after 55 h of incubation at 21° in the presence of increasing IVM doses. Data are presented as mean ± S.D. from 3 to 4 independent experiments. (C) Relative expression of pgp mRNA in IVM-selected C. elegans. mRNA levels were measured by RNA-seq and are expressed as fold change relative to N2 Bristol. Pairwise comparisons between IVM-selected strains and N2 Bristol were performed by using the edgeR package and multiple testing correction was applied with the Benjamini-Hochberg procedure to control the False Discovery Rate (FDR). Data are reported as the mean ± S.D of 3 independent experiments. Genes with a FDR ≤0.05 were considered differentially expressed. ^a^ FDR<0.05; ^b^ FDR<0.01; ^c^ FDR<0.001 vs. wild-type N2 Bristol.Fig. 1

At the end point of the selection process, we compared the sensitivity to IVM of the three resistant strains by larval development assay (Fig. 1B). Compared to the susceptible wild-type N2 Bristol strain, we observed a marked rightward shift in the dose-response curves for IVR10-2014 and IVR10-2022 lines similar to the standard strain IVR10-2009. This fits with a significantly decreased potency of IVM as shown by a more than 6-fold increase in EC_50_ values, i.e., the concentrations at which 50% of the larvae fail to reach adulthood (Table S1). These results show that the three IVM-selected strains acquired resistance with similar dynamics to achieve comparable levels of resistance to IVM.

Divergence of pgp transcriptional expression in IVM-resistant strains

3.2

To assess the impact of IVM selection pressure on the transcriptional profiles of the three IVM-resistant strains we performed RNA-seq analysis and compared the transcriptional expression of pgp genes. Pgp mRNA levels differed relative to N2 Bristol and also among the three resistant lines (Fig. 1C). In the IVR10-2009 strain several pgp were overexpressed, including pgp-5, pgp-11, pgp-13, pgp-6, pgp-7, pgp-4 and pgp-12. The most notable changes concerned pgp-5 (2.8-fold), pgp-6 (4.2-fold), and pgp-7 (7.7-fold). Meanwhile in the IVR10-2014 worms, only pgp-8 and pgp-13 were significantly overexpressed by a 3.9 and 1.7-fold-change, respectively, while pgp-5, pgp-6 and pgp-14 were significantly reduced. Importantly, there was no overlap in the overexpressed pgp genes in the IVR10-2009 and IVR10-2014 strains. Intriguingly, none of the pgp were significantly upregulated in IVR10-2022. Instead, pgp-1 and pgp-6 were significantly downregulated (Fig. 1C).

Given the similar IVM tolerance of the three IVR10 strains, their pgp expression profiles indicate that pgp upregulation is not invariably associated with high IVM tolerance.

Uncoupling of pgp expression and IVM tolerance in IVR10-2009

3.3

Once generated, the IVM-resistant strains were all constantly maintained under IVM selection pressure (10 ng/ml). In order to verify if IVM pressure is required to maintain high tolerance, we removed IVM from the IVR10-2009 NGM plates for 6 months. We compared the sensitivity to IVM measured by larval development assay of IVR10-2009 worms subjected to drug selection pressure (IVR10-2009(+)) with worms maintained on IVM-free NGM plates (IVR10-2009(−)). Both strains showed similar tolerance to IVM (EC_50_ of 9.42 ± 1.31 for IVR10-2009(+) vs. 8.76 ± 2.16 nM for IVR10-2009 (−); Fig. 2A and Table S2). We concluded that long-term maintenance of IVM pressure does not impact IVM tolerance of IVR10 worms.Fig. 2Influence of long-term IVM withdrawal on IVM susceptibility phenotypic resistance and pgp mRNA levels in IVR10–2009.IVM-selected strain IVR10-2009 was maintained for 6 months on selection pressure (IVM at 10 ng/ml) (IVR10-2009(+)), or without IVM (IVR10-2009(−)). (A) Susceptibilities of N2 Bristol and IVM-selected C. elegans strains to IVM assessed by larval development assay. Values represent the percentage of L1 reaching the young adult stage after 55 h of incubation at 21 °C in the presence of increasing doses of IVM. Data are mean ± S.D from 3 independent experiments. **(B)**pgp mRNA levels were measured by RNA-seq. Data are expressed as –fold change relative to N2 Bristol, and reported as the mean ± S.D of 3 independent experiments. Pairwise comparisons between IVM-selected strains and N2 Bristol were performed by using the edgeR package and multiple testing correction was applied with the Benjamini-Hochberg procedure to control the False Discovery Rate (FDR). Data are reported as the mean ± S.D of 3 independent experiments. Genes with FDR ≤0.05 were considered differentially expressed. ^a^ FDR<0.05; ^b^ FDR<0.01; ^c^ FDR<0.001, ^d^ FDR<0.0001 vs. wild-type N2 Bristol.Fig. 2

We then compared pgp transcriptional expression obtained by RNA-seq in worms maintained or not under IVM pressure (IVR10-2009(+) and (−)) (Fig. 2B). Relative to N2 Bristol, IVR10-2009 displayed significant upregulation of several pgp transcripts under IVM selection (+) compared with the unselected condition (−). Compared to previously obtained independent RNA seq data (Fig. 1C), we found an overexpression of pgp-4, pgp-5, pgp-6, pgp-7, pgp-12, and pgp-13. In addition, we observed overexpression of pgp-3, pgp-8, and pgp-10, but no overexpression of pgp-11. The most notable changes concerned pgp-5 (7.7-fold), pgp-6 (7.9-fold), and pgp-7 (11.5-fold). All significant pgp overexpressions observed in the IVR10-2009(+) strain were confirmed by RT-qPCR (Fig. S1). By contrast, only slight pgp overexpression was observed when IVR10-2009 was maintained without IVM pressure (−). This includes overexpression of pgp-10 (2.2-fold), pgp-11 (1.8-fold), and pgp-12 (1.7-fold), whereas the level of expression of pgp-9 was significantly reduced (−1.6-fold). Altogether, these results show that pgp transcriptional upregulation is not necessary to sustain a pre-established resistance to IVM in our experimental conditions. However, keeping the resistant worms under selection pressure could help to maintain a high level of pgp expression, which may have been acquired during the evolutionary selection process.

pgp genes are not induced by acute IVM exposure

3.4

Since IVM resistance was not always associated with constitutive upregulation of pgp genes in our experimental conditions, we tested whether acute IVM exposure could trigger transitory pgp transcriptional overexpression, potentially contributing to the resistance phenotype of the IVR10-2009 strain. We exposed both N2 Bristol and IVR10-2009 C. elegans strains to various concentrations of IVM for 6 h. Following treatment, we assessed the expression levels of pgp genes by RT-qPCR. Relative to each strain's DMSO control, IVM exposure induced no major changes in pgp expression across susceptible and resistant backgrounds (Fig. 3). Indeed, pgp mRNA levels remained unchanged after IVM exposure up to doses of 10 nM and 50 nM for the susceptible and resistant strains, respectively, which correspond to more than 7 × and 5 × their EC_50_ values measured by LDA (Fig. S2). From these concentrations and beyond, only a few pgp genes showed significantly increased expression compared to untreated worms, and the overall level of overexpression remained modest. IVM induced significant upregulation of pgp-10, pgp-11, pgp-12 and pgp-13 in both the N2 Bristol and IVR10-2009 strains, and of pgp-14 in IVR10-2009 with fold-changes ranging from 2 to 3 only. Furthermore, expression of pgp does not exhibit a dose-response profile to IVM. For example, pgp-10 was significantly upregulated in the N2 Bristol at 10 ng/ml of IVM but its expression decreased at higher doses. Finally, pgp-13 exhibited a significant downregulation, both in IVR10-2009 and N2 Bristol at the highest IVM doses (Fig. 3, Fig. S2).Fig. 3Impact of IVM acute exposure on pgp transcription profile.L4 worms were exposed 6 h to 2.5 to 30 nM IVM for the N2 Bristol (A) and 7.5 to 100 nM IVM for the IVR10-2009 (B). RNAs were extracted and levels were measured by real-time qPCR. Data are expressed as –fold change relative to the DMSO condition, and reported as the mean of 3 independent experiments. Changes in mRNA levels were normalized with respect to *tba-*1 mRNA levels. Statistical significance between the control (DMSO) and IVM exposed conditions are shown in Fig. S2.Fig. 3. Fig. 4Temporal transcription profile of pgp expression following withdrawal of IVM pressure.****(A) The IVR10-2009 strain was maintained for up to 6 months either under IVM pressure (10 ng/mL IVM; denoted IVR10(+)) or without selection (0 ng/mL IVM; denoted IVR10(−)). Following this initial phase, defined as time point M0, each condition was either continued or switched to the opposite treatment (±IVM) for an additional 6 months, resulting in four selection pressure histories: IVR10(+,+); IVR10(+,-); IVR10(-,+) and IVR10(-,-). pgp mRNA levels were measured by real-time quantitative PCR (RT-qPCR) at M0, M3 (after 3 months), and M6 (after 6 months). (B) Transcriptional expression data are presented as fold change relative to N2 Bristol. Statistical significances compared to the M0 condition are shown in Fig. S4. Data are reported as the mean of 3 independent experiments.Fig. 4

All together, these results suggest that ML resistance in C. elegans does not result from pgp induction after acute drug exposure.

pgp overexpression is transient and not induced by prolonged IVM exposure in the IVR10-2009 worms

3.5

We next took advantage of our previous observation that pgp transcriptional overexpression was decreased and returned near N2 Bristol levels after withdrawing IVM pressure for 6 months (Fig. 2.) to precisely investigate the dynamics of pgp expression in relation to IVM exposure. We compared pgp expression in IVM-selected IVR10-2009(+) and IVR10-2009(−) strains in response to prolonged IVM exposure or withdrawal, over a long timescale, ranging from weeks to several months. To this end, the IVR10-2009(+) strain was maintained for up to six months either with (IVR10-2009(+,+)) or without IVM pressure (IVR10-2009(+,-)) (Fig. 4A). In parallel, IVR10-2009(−) line previously cultured without IVM for 6 months, was further maintained continuously without selection pressure (IVR10-2009(-,-)), or was re-exposed to IVM (IVR10-2009(-,+)) (Fig. 4A). Then, we monitored pgp transcript levels by RT-qPCR every week for a month after the switch (Fig. S3). No major pgp deregulation emerged, indicating that any potential changes may occur progressively over a longer period of time. We then assessed whether IVM exposure could influence pgp expression over a longer period of up to 6 months. First, after IVM removal for 6 months, we confirmed that the IVR10-2009(+) strain failed to maintain pgp genes overexpression as the expression level of pgp-5, pgp-6 and pgp-7 decreased significantly by a factor of 2 to 6 (Fig. 4B and Fig. S4). By contrast, when IVM pressure was maintained, only pgp-7 showed significantly decreased expression. Importantly, in the IVR10-2009(−) line, which no longer exhibited pgp overexpression, we observed no significant changes in pgp mRNA levels after re-exposure of worms to IVM (Fig. 4B, Fig. S4). Taken together, these results suggest that prolonged IVM exposure does not result in induction of pgp genes but may allow some level of overexpression to be maintained.

pgp upregulation is associated with early steps of IVM resistance emergence in C. elegans

3.6

To better understand how pgp gene transcriptional expression evolves during the progressive acquisition of IVM resistance, we monitored transcript levels at several key stages of the selection process in the IVR10-2022 strain (Fig. 5A) and in the IVR10-2014 strain (Fig. 5B). pgp expression was quantified by RT-qPCR at representative time points corresponding to concentration plateaus: early (2a) and late (2b) stages at 2 ng/mL, an intermediate step at 6 ng/mL (6), and the final plateau at 10 (10) or 11 ng/mL (11).Fig. 5pgp expression profile during IVM selection of IVR10–2022 and IVR10–2014 resistant C. elegans strains.Selection curves representing the dynamic process leading to the generation of IVR10-2014 (A) and IVR10-2022 (B) strains, as in Fig. 1A. The time points at which pgp gene transcriptional expression was analyzed are indicated and correspond to specific stages of the selection process: early stage of the 2 ng/mL plateau (2a), end of the 2 ng/mL plateau (2b), 6 ng/mL (6), and the end of the last plateau at 10 ng/ml (10) for IVR10-2014 and 11 ng/mL (11) for IVR10-2022. These samples represent key transitional phases in the acquisition of resistance under increasing IVM pressure. Dynamics of the pgp transcript levels during the selection process was measured by real-time quantitative PCR (RT-qPCR) for IVR10-2014 (C) and IVR10-2022 (D). Expression data are presented as fold change relative to the N2 Bristol at the start of the selection, which was normalized to 1. Values represent the mean of 3 independent experiments. Transcript levels were normalized to *cdc-*42 mRNA expression. Statistical significance between the control (DMSO) and IVM exposed conditions is shown in Fig. S5.Fig. 5

Focusing on the early plateau at 2 ng/mL, both lines already showed significant pgp induction at stage 2a: IVR10-2014 for pgp-5, pgp-6, pgp-7, and pgp-8 (3 to 5-fold), and IVR10-2022 for pgp-5, pgp-6, pgp-9, pgp-10, and pgp-14 (3 to 8-fold). At stage 2b, IVR10-2014 displayed a decline in pgp overexpression with only modest upregulation for pgp-5, pgp-6, and pgp-14 (1.5 to 2-fold). In contrast, IVR10-2022 showed a dramatic peak at stage 2b for pgp-5 (12-fold) and pgp-6 (13-fold), as well as sustained upregulation for other pgp already overexpressed at stage 2a (Fig. 5C–D, Fig. S5). Beyond the initial plateau at 2 ng/ml, pgp overexpression showed a progressive attenuation through later stages even though the level of tolerance to IVM increased sharply, underscoring that the most pronounced transcriptional shifts occur at the 2 ng/mL plateau (Fig. 5C–D, Fig. S5).

Together, these data indicate that progressive adaptation to IVM is predominantly accompanied by strong early pgp induction at the 2 ng/mL plateau, supporting a role for PGPs in the initial phase of resistance acquisition.

The transcriptomic program of resistance is enriched in ciliary and neuronal pathways

3.7

In order to assess the variation in gene expression in resistant lines beyond pgp genes, we performed RNA-seq analyses in the resistant IVR10-2009, IVR10-2014, IVR10-2022 and the wild-type N2 Bristol strains. Compared with the N2 Bristol strain, the resistant lines exhibited a variable number of differentially expressed genes (DEGs) with a FDR <0.05 (Table S3). The IVR10-2022 strain exhibited only 671 DEGs while we found 5484 and 7264 DEGs in the IVR10-2009 and the IVR10-2014 strains respectively (Fig. 6A). This high variability in gene expression of worms that acquired an IVM resistance phenotype could reflect the evolutionary history of stochastic resistance acquisition. Nevertheless, we observed that the three IVM resistant strains had 359 DEGs in common, suggesting that they might share similar resistance acquisition mechanisms, at least partly. To gain insight into the functional impact of gene expression changes associated with IVM resistance, we performed enrichment analyses on the differentially expressed genes shared by all the resistant C. elegans strains. We assessed enrichment by phenotypic annotations (Fig. 6B), biological processes (Fig. 6C), and tissue-specific expression (Fig. 6D). Interestingly, resistance-associated DEGs were enriched in biological processes related to intraciliary transport and cilium organization, many of which are known to confer the Dyf phenotype and IVM resistance (Brinzer et al., 2024; Page, 2018), including bbs-8, che-3, che-11, dyf-1, dyf-2, dyf-6, dyf-11, dyf-13, osm-6 and osm-12 (Table S4). Unsurprisingly, dye filling defect was one of the most significantly enriched phenotype along with IVM or anthelmintic response variant and IVM resistance. Moreover, enrichment of DEGs within ciliated, phasmid and amphid neurons tissues was also prominent.Fig. 6**Functional enrichment analysis of differentially expressed genes in various resistant C. elegans strains: biological processes, phenotypes, and anatomical expression.**Differentially expressed genes (DEGs) of 3 independent experiments were identified from RNA-seq datasets generated for each resistant strain and correspond to genes significantly (FDR <0.05) up- or downregulated compared to the susceptible strain N2 Bristol. Gene set enrichment analyses were then performed. (A) The Venn diagram illustrates the overlap of differentially expressed genes between the resistant strains. Lollipop plots were generated by using ShinyGo web application and display significantly enriched terms in (B) phenotypes, (C) gene ontology biological processes and (D) tissue-specific expression domains.Fig. 6

Thus, despite divergent transcriptomic profiles reflecting distinct evolutionary histories of resistance acquisition, the enrichment analyses consistently point towards shared modulation of ciliary amphid function associated with IVM resistance. This highlights amphid sensory structures as a central hub in the molecular response to IVM pressure.

The development of IVM resistance is accompanied by the acquisition of the Dyf phenotype

3.8

The red, lipophilic tracer DiI is used to label ciliated sensory neurons and to assess structural or functional defects. Both IVR10-2009 and IVR10-2014 have been reported as Dyf, an established marker of amphid dysfunction. Here, we show that IVR10-2022 is also Dyf (Fig. S6), supporting the view that impaired ciliated sensory neuron function may be a prerequisite for ML resistance.

To better understand the relationship between the Dyf phenotype and IVM resistance, we monitored the integrity of ciliated sensory neurons by Dil staining at different stages of the selection process in IVR10-2014 and IVR10-2022 strains (Fig. 5A and B) (Fig. 7A and B). Early in selection, amphid neurons readily absorbed DiI. However, labeling was completely lost from the end of the 2 ng/ml plateau in IVR10-2014 or, at least, from the 6 ng/ml dose in IVR10-2022. Interestingly, these steps in the selection process coincide with a sharp increase in resistance to high IVM doses, suggesting that the Dyf phenotype is tightly linked to the establishment of ML resistance.Fig. 7. Amphid dye uptake.Adult worms from the selection process of IVR10-2014 (A) or IVR10-2022 (B) and corresponding to time points: early stage of the 2 ng/mL plateau (2a), end of the 2 ng/mL plateau (2b), 6 ng/mL (6), and the end of the last plateau at 10 ng/ml (10) for IVR10-2014 and 11 ng/mL (11) for IVR10-2022 (Fig. 5A and B) were incubated with the red fluorescent lipophilic dye Dil.Fig. 7

Discussion

4

PGPs have long been implicated in modulating IVM susceptibility, and their transcriptional overexpression is frequently reported in ML–resistant parasites. However, their precise contribution to the resistant phenotype has remained unsolved. In our study, the transcriptomic profiling of multiple ML-resistant C. elegans strains revealed highly heterogeneous pgp expression patterns, arguing against a consistent correlation between pgp upregulation and stable IVM resistance. Our data instead suggest that transient pgp induction may contribute to early tolerance, facilitating the subsequent selection of more robust resistance mechanisms. In line with this, enrichment analyses highlighted a strong and recurrent deregulation of neuronal and ciliary gene networks across resistant strains. These signatures converge with the Dyf phenotype, which develops for worms able to grow in the presence of IVM concentrations lethal to susceptible strains.

A key strength of our study lies in the use of three distinct IVM-resistant strains: IVR10-2009, IVR10-2014, and IVR10-2022, each developed independently by stepwise selection on increasing concentrations of IVM. This experimental-evolution approach provides a more realistic view of how IVM resistance develops than mutagenesis-based models. After ∼80–90 generations, all strains converged on equivalent drug-resistance levels, following comparable dynamics. All selection processes acquired linear IVM tolerance interspersed with selection plateaus, during which populations adapted before tolerating higher doses. A 2 ng/mL plateau was common to all strains but varied in duration, notably lasting ∼20 weeks in IVR10-2022. A second plateau at 6–7 ng/mL occurred in IVR10-2009 and IVR10-2014 but was absent in IVR10-2022. These patterns suggest that independent C. elegans populations adapt via shared mechanisms with lineage-specific contingencies, shaped by exposure history, stochastic events, or genetic background (Baran et al., 2024; Estes et al., 2011; Hellinga et al., 2024; McGrath et al., 2011). Here, the use of three distinct IVM-resistant strains allowed us to identify both shared traits (e.g., dye-filling defects, resistance to 10–11 ng/mL IVM) and strain-specific features, such as pgp transcriptional expression profiles.

Our transcriptomic study performed on the IVM-selected strains at the end-point of the selection process showed that pgp expression differed across the three lineages. We observed a marked upregulation of several pgp genes confined to only one resistant strain (IVR10-2009), while discrete overexpression occurred in 2014, and none in 2022. This heterogeneity is consistent with prior findings in C. elegans and parasitic nematodes, where pgp expression under ML pressure is context-dependent, and often lineage-specific (Lespine et al., 2024). Indeed, not all resistant isolates show pgp upregulation, and in some cases, specific pgp genes are even downregulated (Lespine et al., 2024). These findings may reflect both the plasticity of adaptive responses and the functional redundancy among pgp paralogs (Langeland et al., 2021; Lespine et al., 2024), many of which are differentially expressed across tissues (Ghaddar et al., 2023), and they further call into question the extent to which PGPs are required to sustain IVM resistance once it has been established.

The IVR10-2009 lineage has been previously characterized by James and Davey (2009) (James and Davey, 2009) and Ménez et al. (2016) (Ménez et al., 2016), who showed overexpression of pgp-1 (both studies), and pgp-6, pgp-9, and pgp-14 (Ménez et al. only). In our study, this lineage showed both conserved and divergent regulation: pgp-6 remained upregulated, but pgp-1 and pgp-14 expressions were comparable to N2 Bristol. These observations suggest that the level of pgp expression could fluctuate during maintenance of resistant worms in culture. This notion is further reinforced by the IVR10-2009(−) line, which, despite being maintained without continuous IVM pressure, retained a stable resistant phenotype while showing little to no pgp overexpression. In addition, these findings must be interpreted cautiously, as PGPs function depends not only on transcription but also on protein stability, folding, trafficking, and ATPase activity (Buccitelli and Selbach, 2020; Kramer et al., 1995; Liu et al., 2014; Srivastava et al., 2022; Zhang et al., 2004). Complementary assays such as inhibitor sensitivity, fluorescent drug efflux assays, or tissue-specific pgp expression, are needed to confirm PGP activity. Nonetheless, pgp genes have been implicated in IVM tolerance, as pgp knockout mutants show hypersensitivity in both wild-type and IVM-resistant backgrounds (Ardelli and Prichard, 2013; Blancfuney et al., 2025; Janssen et al., 2013; Ménez et al., 2016), while pgp overexpression has been shown to protect C. elegans against IVM (Gerhard et al., 2021; Janssen et al., 2015; Ménez et al., 2016).

In our dataset, pgp-5, pgp-6, and pgp-7 were overexpressed in only one strain (IVR10-2009), while pgp-8 was overexpressed in another (IVR10-2014). Notably, all these genes are located on chromosome X. The strain-specific overexpression of distinct pgp genes suggests that PGP-mediated drug tolerance is flexible and not conserved. However, their shared chromosomal location may indicate common regulatory control, but their differential activation suggests context-dependent usage. According to prior single-cell RNA-seq data, pgp-5 and pgp-6 are expressed in the intestine and neurons, pgp-7 in the tail (weakly), and pgp-8 in phasmid neurons (Ghaddar et al., 2023). These tissues represent plausible entry or clearance sites for xenobiotics, supporting a model in which drug tolerance may involve tissue-specific efflux rather than systemic detoxification. What remained unclear was whether pgp transcriptional expression is induced in response to acute IVM exposure and once the resistance is stabilized.

Several studies in C. elegans have shown conflicting results on the effect of IVM exposure on pgp transcriptional expression with studies showing no induction upon IVM exposure while others reported upregulation of some pgp (Ardelli and Prichard, 2013; Bygarski et al., 2014; Dube et al., 2023; Laing et al., 2012). To address these discrepancies, we performed short-term IVM exposure both on the N2 Bristol and IVR10-2009 strains. We observed only moderate upregulation of pgp-10 to pgp-14, with no effect on pgp-1 to pgp-9, for doses of IVM well above the EC_50_ for each of the lines. This agrees with previous studies showing limited pgp induction even at IVM doses 100-fold above the EC_50_ for N2 Bristol (Dube et al., 2023) or in the triple GluCl mutant strain (Laing et al., 2012), suggesting that pgp gene regulation is not part of an immediate stress response to IVM.

To address the role of pgp transcriptional expression in long-term resistance, we examined transcriptional and phenotypic stability of the IVR10-2009 strain under sustained IVM exposure or following withdrawal. We observed that pgp transcript levels declined over time in the absence of the drug, while resistance remained stable. This finding mirrors prior results in C. elegans, where pgp expression dropped after IVM removal (James and Davey, 2009). Re-exposure to IVM failed to restore pgp expression, supporting the idea that pgp induction is transient, non-reinducible, and mainly linked to early selection. We addressed this latter point by monitoring pgp expression throughout the selection process. Our data from the IVR10-2014 and IVR10-2022 strains provide direct evidence that early tolerance is associated with a common upregulation of pgp genes, notably pgp-5 and pgp-6. In addition, strain-specific responses were observed: pgp-7 and pgp-8 were induced exclusively in IVR10-2014, whereas pgp-9, pgp-10, and pgp-14 peaked specifically in IVR10-2022 at the 2 ng/mL plateau. These observations support a protective role of PGPs against the harmful effects of IVM during the very first stages of selection. In C. elegans, overexpression of pgp-6 increases tolerance to IVM (Ménez et al., 2019), while expression of Parascaris sp. pgp-11 driven by the C. elegans pgp-11 promoter reduces sensitivity to IVM (Janssen et al., 2015). Similarly, intestinal expression of P. univalens pgp-9 confers a protective effect against the drug (Gerhard et al., 2021). Therefore, our observations are in agreement with the idea that pgp activation could facilitate early adaptation to moderate IVM exposure. This hypothesis is consistent with the concept that repeated exposure to subtherapeutic doses can promote the selection of resistance. Conversely, if such doses do not permit survival, selection of resistance cannot occur (Coles et al., 2005; James and Davey, 2009; Ménez et al., 2016, 2019; Van Zeveren et al., 2007). Accordingly, nhr-8 mutants, which are hypersensitive to IVM, fail to develop resistance under stepwise exposure. This likely reflects the role of NHR-8 in regulating xenobiotic detoxification, including pgp expression (Ménez et al., 2019). While we did not examine nhr-8 here, the pgp expression pattern we report here is compatible with sustained pgp expression early in the selection process to tolerate moderate doses of IVM and facilitate resistance development. Besides, pgp upregulation did not persist in the IVR10-2014 and IVR10-2022 strains despite a dramatic increase in resistance suggesting that pgp activation reflects a costly, early-phase adaptation later replaced by more stable mechanisms. Altogether, our transcriptional data indicate that pgp overexpression could promote resistance emergence but is dispensable for long-term maintenance. Distinguishing early from stable mechanisms is essential for interpreting anthelmintic resistance trajectories.

We focused primarily on pgp genes to investigate transcriptomic regulation during the early stages of the selection process. While many studies have highlighted the role of pgp genes in ML tolerance through drug efflux, additional and previously unsuspected factors may also contribute to worm survival at moderate IVM doses. For example, members of the nspc gene family show extensive transcriptional changes, with 13 out of 20 genes significantly dysregulated in resistant lines compared to the N2 Bristol strain (Table S3). However, the functional roles of nspc genes remain largely unknown, preventing us from determining whether they could contribute to IVM resistance. A more comprehensive transcriptomic analysis of lines at intermediate selection stages, combined with improved characterization of poorly studied genes, may therefore uncover additional mechanisms underlying adaptation to IVM and open new research avenues. In addition, we cannot exclude the possibility that epigenetic processes contribute to early adaptation. Several studies have shown that environmental stresses can induce transient, transgenerational epigenetic modifications (Filipowicz and Allard, 2025), which could potentially underlie the transcriptional changes we observed for pgp genes during the initial phases of selection.

To identify the mechanisms underlying stable resistance to ML, we conducted an RNA-seq study on the three resistant strains. Although the number of deregulated genes was extremely variable between resistant lines, we found 359 genes differentially expressed compared to N2 Bristol that were common to all strains. Gene set enrichment analyses identified consistent deregulation of ciliary gene networks in sensory neurons such as amphids across resistant strains, pointing to a shared sensory-based adaptation to drug pressure. These transcriptomic signatures are consistent with the Dyf phenotype, an established marker of amphid dysfunction, previously reported to confer IVM resistance in C. elegans (Page, 2018). Interestingly, both the IVR10-2009 and 2014 strains exhibit this phenotype (James and Davey, 2009; Ménez et al., 2016). In addition, we showed that the IVR10-2022 lineage is also Dyf, suggesting that the Dyf phenotype is tightly linked to the establishment of ML resistance. Interestingly, worms became Dyf at the end of the 2 ng/ml IVM plateau for strains IVR10-2014 and 2022, just as resistance was increasing sharply, supporting that resistance establishment is closely associated with loss of amphid functionality. Amphid defect may also play a role in resistance to ML in nematode parasites as amphid structural changes and Dyf phenotype have been described in ML-resistant H. contortus isolates, with evidence of shortened or degraded amphids (Freeman et al., 2003; Urdaneta-Marquez et al., 2014). In C. elegans, most of the genes that confer both the Dyf phenotype and IVM resistance carry nonsense mutations (Brinzer et al., 2024; Page, 2018). This is exemplified by the IVR10-2009 line, in which a mutation in dyf-7 largely accounts for both traits (Urdaneta-Marquez et al., 2014). Our transcriptomic analysis further shows that Dyf and IVM-resistance genes are also strongly deregulated at the transcriptional level, indicating that selection under IVM imposes both mutational and transcriptional pressures. Further investigations will be necessary to assess whether these transcriptional changes contribute directly to IVM resistance, or instead reflect compensatory or reorganizational responses to amphid structural defects caused by genetic mutations remains to be determined.

The Dyf phenotype alone may not always be sufficient to confer strong ML resistance, as the IVR10-2014 line still required 10 weeks to withstand IVM doses above 7 ng/ml, despite already carrying this trait. Still, many studies show that mutations conferring the Dyf phenotype led to IVM-resistance but how is this achieved remains to be determined. It has been proposed that altered amphidial architecture may lead to decreased drug uptake (Page, 2018) but the picture seems to be more complex several lines of evidence indicate that drug resistance in Dyf mutants cannot be solely explained by altered permeability to MLs. First, pgp overexpression in intestinal and cuticular tissues reduces IVM sensitivity (Gerhard et al., 2021), suggesting multiple routes of ML entry. Additionally, Guerrero et al. (2021) observed upregulation of pgp genes in a Dyf mutant, proposing that amphidial defects and efflux mechanisms may act synergistically. Lately, Li et al. showed that mutants of ubr-1, an E3 ubiquitin ligase, which exhibit Dyf phenotype, mediate aberrant glutamate signaling leading to IVM resistance in C. elegans trough IVM-targeted GluCls downregulation. Further studies are therefore necessary to precisely dissect the mechanisms linking the functionality of amphids to resistance to MLs.

Conclusion

5

Our findings challenge the hypothesis that pgp transcriptional overexpression is a sustained core mechanism of IVM resistance in stabilized resistant populations. First, resistance can persist without maintaining pgp upregulation. Second, pgp are not inducible by acute or chronic IVM exposure in our C. elegans model of IVM resistance. Rather, pgp transcriptional expression peaks during early adaptation under moderate drug pressure. These data support a model in which early pgp activation aids survival during selection, facilitating emergence of more stable resistance mechanisms. Among these, alterations in amphidial structures and ciliary pathways are recurrent features across resistant strains. We propose a two-step model: a transient PGP-based tolerance phase under moderate IVM pressure, followed by stable resistance based on sensory neuron remodeling. Understanding how the Dyf phenotype accompanies the resistant state will be important to illuminate the upstream molecular mechanisms that give rise to this adaptation.

CRediT authorship contribution statement

Lucas Barat: Writing – review & editing, Writing – original draft, Visualization, Validation, Investigation, Formal analysis, Data curation, Conceptualization. Eva Guchen: Writing – review & editing, Investigation. Clara Blancfuney: Writing – review & editing, Investigation. Marie Garcia: Writing – review & editing, Investigation. Mélanie Alberich: Writing – review & editing, Investigation. Anne Lespine: Writing – review & editing, Project administration, Funding acquisition. Rémy Betous: Writing – review & editing, Writing – original draft, Visualization, Validation, Supervision.

Funding

This work was funded by the French National Research Agency (ANR, DePAR Contract n°ANR-22-CE44-0022), and France Futur Elevage (F2E) Institut-Carnot Santé Animale (ANTHERIN Contract n° 00005255).

Conflict of interest

The authors declare no conflict of interest.