# Spatio-temporal diversity and genetic architecture of pyrantel resistance in Cylicocyclus nassatus, the most abundant horse parasite

Guillaume Sallé, Élise Courtot, Cédric Cabau, Hugues Parrinello, Delphine Serreau, Fabrice Reigner, Amandine Gesbert, Lauriane Jacquinot, Océane Lenhof, Annabelle Aimé, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Guillaume Sallé, Élise Courtot, Cédric Cabau, Hugues Parrinello, Delphine Serreau, et al. Spatio-temporal diversity and genetic architecture of pyrantel resistance in Cylicocyclus nassatus, the most abundant horse parasite. Peer Community Journal, 2025, 5, pp.e69. 10.24072/pcjournal.571. hal-05197483

HAL Id: hal-05197483

https://hal.inrae.fr/hal-05197483v1

Submitted on 2 Aug 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Research article

Published
2025-07-07

**Cite as**
Guillaume Sallé, Élise Courtot, Cédric Cabau, Hugues Parrinello, Delphine Serreau, Fabrice Reigner, Amandine Gesbert, Lauriane Jacquinot, Océane Lenhof, Annabelle Aimé, Valérie Picandet, Tetiana Kuzmina, Oleksandr Holovachov, Jennifer Bellaw, Martin K. Nielsen, Georg von Samson-Himmelstjerna, Sophie Valière, Marie Gislard, Jérôme Lluch, Claire Kuchly and Christophe Klopp (2025) *Spatio-temporal diversity and genetic architecture of pyrantel resistance in Cylicocyclus nassatus, the most abundant horse parasite*, Peer Community Journal, 5: e69.

gsalle.inra.email@gmail.com

**Peer-review**
Peer reviewed and
recommended by
PCI Genomics,
https://doi.org/10.24072/pci-
genomics.100290

<figure><img src="image_3.png" /></figure>

BY

This article is licensed
under the Creative Commons
Attribution 4.0 License.

# Spatio-temporal diversity and genetic architecture of pyrantel resistance in *Cylicocyclus nassatus*, the most abundant horse parasite

Guillaume Sallé<sup>1,1</sup>, Élise Courtot<sup>1,1</sup>, Cédric Cabau<sup>2</sup>, Hugues Parrinello<sup>1,3</sup>, Delphine Serreau<sup>1</sup>, Fabrice Reigner<sup>1,4</sup>, Amandine Gesbert<sup>4</sup>, Lauriane Jacquinot<sup>1</sup>, Océane Lenhof<sup>1</sup>, Annabelle Aimé<sup>5</sup>, Valérie Picandet<sup>5</sup>, Tetiana Kuzmina<sup>1,6,7</sup>, Oleksandr Holovachov<sup>1,8</sup>, Jennifer Bellaw<sup>9</sup>, Martin K. Nielsen<sup>9</sup>, Georg von Samson-Himmelstjerna<sup>1,10</sup>, Sophie Valière<sup>11</sup>, Marie Gislard<sup>11</sup>, Jérôme Lluch<sup>11</sup>, Claire Kuchly<sup>1,11</sup>, and Christophe Klopp<sup>12</sup>

Volume 5 (2025), article e69

https://doi.org/10.24072/pcjournal.571

## Abstract

Cyathostomins are a complex of 50 intestinal parasite species infecting horses and wild equids. The massive administration of modern anthelmintic drugs has increased their relative abundance in horse helminth communities and selected drug-resistant isolates worldwide. *Cylicocyclus nassatus* is the most prevalent and the most abundant species. The tedious identification and isolation of these worms have hampered studies of their biology that remain largely uncharacterised. Here we have leveraged ultra-low input sequencing protocols to build a reference genome for the most prevalent horse strongyle species. Using this resource, we have established the first estimates of its genetic diversity and population structure on a gradient ranging from Ukraine (close to modern horse domestication area) to North America, while capturing a 19th-century snapshot of *C. nassatus* diversity in Egypt. Our results support a diverse and lowly structured global population. Modern populations displayed lower nucleotide diversity relative to the old North African isolate. We identified the first genetic candidates upon which pyrantel (an anthelmintic drug used in companion animals) selection likely applied in field populations, highlighting previously suspected genes coding for nicotinic acetylcholine receptor subunits, and identifying new candidates showing differential expression in independently evolved *Caenorhabditis elegans* lines. These results offer a first resource to widen current knowledge on cyathostomin biology, unravel novel aspects of pyrantel resistance mechanisms and provide candidate genes to track pyrantel resistance in the field.

<footnote n="1"><sup>1</sup>Université de Tours, INRAE UMR1282 Infectiologie et Santé Publique, Nouzilly - France, <sup>2</sup>Sigenae, GenPhySE, Université de Toulouse, INRAE, ENVT, Castanet Tolosan - France, <sup>3</sup>MGX-Montpellier GenomiX, Univ. Montpellier, CNRS, INSERM, Montpellier - France, <sup>4</sup>INRAE UE 1297 Unité Expérimentale de Physiologie Animale de l'Orfrasière, Nouzilly - France, <sup>5</sup>Centre hospitalier Vétérinaire Équin de Livet, 1497 Route de Castillon, Livarot Pays d'Auge- France, <sup>6</sup>I. I. Schmalhausen Institute of Zoology NAS of Ukraine, Kyiv - Ukraine, <sup>7</sup>Institute of Parasitology, Slovak Academy of Sciences, Košice - Slovak Republic, <sup>8</sup>Department of Zoology, Swedish Museum of Natural History, Stockholm - Sweden, <sup>9</sup>Department of Animal and Veterinary Sciences, Aarhus University, Tjele, Denmark, <sup>10</sup>Institute for Parasitology and Tropical Veterinary Medicine, Freie Universitaet Berlin, Robert-von-Ostertag-Str. 7, Berlin - Germany, <sup>11</sup>INRAE, US 1426, GeT-PlaGe, Genotoul, France Génomique, Université Fédérale de Toulouse, Castanet Tolosan - France, <sup>12</sup>Sigenae, Genotoul Bioinfo, MIAT UR875, INRAE, Castanet Tolosan - France</footnote>
Introduction

Parasitic nematodes account for 11.8 million disability-adjusted life years in humans (DALYs &
Hale Collaborators, 2017) and inflict major losses on pet and livestock species (Kaplan &
Vidyashankar, 2012; von Samson-Himmelstjerna et al., 2021; Nielsen, 2022; Marsh & Lakritz, 2023).
While prevention and sanitation have played a major role in bringing the Guinea worm to the brink of
extinction (Durrant et al., 2020), parasite control programs still heavily rely on the administration of
anthelmintic drugs in both the medical (Bradley et al., 2021; Gandasegui et al., 2022) and veterinary
settings (Laing et al., 2017; von Samson-Himmelstjerna et al., 2021; Nielsen, 2022; Marsh & Lakritz,
2023). In livestock-infecting species, anthelmintic treatments have reduced nematode diversity at both
the species and community levels and species scales as illustrated in horses. In this rich host-parasite
system, the once predominant horse large strongyles (Strongylus spp.) are now encountered at low
prevalence across managed horse populations (Jürgenschellert et al., 2022) while they dominate
parasite communities in a Canadian feral horse population (Jenkins et al., 2020). Field observations
also support increased odds to encounter S. vulgaris on farms implementing evidence-based
treatments (Tydén et al., 2019). On the contrary, cyathostomins, a complex of fifty described species
inhabiting the horse hindgut and dominated by Cylicocyclus nassatus (Ogbourne, 1972; Lyons et al.,
1992; Bucknell et al., 1995; Collobert-Laugier et al., 2002; Kuzmina et al., 2016), have become
predominant (Herd, 1990) and represent the most important cause of parasite-mediated death in
young horses (Sallé et al., 2020). Cyathostomins are encountered worldwide across a wide range of
equids, including donkeys, horses and zebras (Kuzmina & Kuzmin, 2008; Lichtenfels et al., 2008;
Kuzmina et al., 2013; Tombak et al., 2021). However, investigation of their biology has been hampered
by the complexity of their assemblages that usually encompass more than ten species within a horse,
and the tediousness of their morphological identification (Lichtenfels et al., 2008; Bellaw & Nielsen,
2020). Recent applications of metabarcoding experiments (Poissant et al., 2021) have uncovered
additional aspects of their response to drugs (Nielsen et al., 2022) and plant products (Malsa et al.,
2022), their interaction with their host gut microbiota (Boisseau et al., 2023) or fluctuation in their
relative abundance throughout a grazing season in unmanaged horses (Sargison et al., 2022). But the
genetics of drug resistance in cyathostomins remains unknown for drugs other than benzimidazoles
(Hodgkinson et al., 2008).

Beyond reshaping helminth community structure, anthelmintic drugs have also remodelled the
genetic diversity of parasite species as evidenced in *Haemonchus contortus* (Sallé et al., 2019; Doyle
et al., 2019, 2022), a blood feeding parasite of small-ruminants. In that case, drastic loss of genetic
diversity occurred over a beta-tubulin coding gene in benzimidazole-resistant isolates (Sallé et al.,
2019). Additional efforts of breeding experimental back-cross lines have also shed light on the
molecular architecture of resistance to levamisole (Doyle et al., 2022a), ultimately yielding genetic
markers to anticipate the selection of resistant isolates in the field (Antonopoulos et al., 2022).
Altogether, the role of candidate genes already proposed to affect anthelmintic sensitivity to
benzimidazoles (Kwa et al., 1995) or levamisole (Neveu et al., 2010) have been confirmed using
genome-wide mapping approaches, but discrepancies have emerged for the genetic structure of
ivermectin resistance (Laing et al., 2022). In addition, the genes associated with the resistance to other
drugs like pyrantel remain poorly characterised.

Pyrantel is a widely used anthelmintic drug for parasite control in humans (Moser et al., 2017), pets (Kopp et al., 2008) and horses (Lyons et al., 1999). Reduced efficacy has been observed in the dog hookworm *Ancylostoma caninum* and frequently occurs in the horse cyathostomins (Lyons, 2003; Sallé et al., 2017). Early electrophysiology measures in muscle vesicles of the pig parasites *Oesophagostomum dentatum* (Robertson et al., 1994) or *Ascaris suum* (Robertson et al., 2000) concluded that pyrantel acts as an acetylcholine agonist on cation channels present at neuromuscular
junctions. This ultimately results in their spastic paralysis and elimination from their host. The knock-
down of the genes coding for some of these receptor subunits was affecting *C. elegans* susceptibility
towards pyrantel, either completely for *unc-29* and *unc-63* knock-outs or partially for *unc-38* knock-
outs (Sleigh, 2010). In addition, pyrantel-resistant *A. caninum* isolated from dogs had reduced
expression of these genes in comparison to a susceptible isolate (Kopp et al., 2009). Other
heterologous expression experiments of *H. contortus acr-8* (Blanchard et al., 2018) or *acr-26* and *acr-
27* (*H. contortus* and *Parascaris* sp.) in *C. elegans* also increased the resistance level to this drug
(Courtot et al., 2015). Despite these converging strands of evidence, it is yet unknown whether pyrantel
treatment in the field primarily select variants of these candidates, whether additional candidates
contribute to decreased pyrantel sensitivity and if the known targets shared by *C. elegans* and *A.
caninum* are under selection in cyathostomin populations.

To bridge this knowledge gap, this study aims to investigate the genetic architecture of pyrantel
resistance in cyathostomin populations using a set of pyrantel-sensitive and -resistant isolates. The
success of mapping approaches to isolate drug resistance candidates has relied on well-annotated
chromosome-level assemblies. While a substantial number of assemblies have recently been
produced for helminth species of medical and veterinary importance (International Helminth Genomes,
2019), genome assembly has proved challenging for parasitic nematodes due to the limited genetic
material available for each individual and their highly heterozygous and repetitive medium-sized
genomes (International Helminth Genomes, 2019; Doyle, 2022). In some cases like the horse
cyathostomins, the complexity to collect, isolate and identify worm material from horses (Lichtenfels
et al., 2008) adds additional challenges to produce sufficient quantities of good quality DNA and RNA
(Louro et al., 2021). To date, *Cylicostephanus goldi* is the only cyathostomin species with a released
genome (International Helminth Genomes, 2019) and transcriptome (Cwiklinski et al., 2013)
assemblies. However, short-read data generated for this species could not resolve the full genome
sequence and yielded a heavily fragmented assembly (International Helminth Genomes, 2019).
Nonetheless, the recent description of ultra-low input protocol applicable to long-read sequencing
technologies (Kingan et al., 2019) opens new perspectives for the generation of high-quality genome
assemblies and subsequent study of genetic diversity for hard-to-collect species.

Using this approach and Hi-C genomic analysis technique, we have built a chromosomal assembly
from a single *Cylicocyclus nassatus* individual and applied Hi-C technology to scaffold a chromosome-
level genome. We investigated the diversity of this species using contemporary worms covering an
East-West gradient from Ukraine to Kentucky, USA and established patterns of variation between
these contemporary populations and an old isolate collected in Egypt in the 19th century. Finally,
comparison of the genetic diversity of pyrantel-resistant and -susceptible isolates identified novel
candidates for pyrantel resistance in *C. nassatus*. Orthologs of these candidate genes were found
differentially expressed in evolved *C. elegans* lines.

**Material and methods**

Worm collection procedures

Unless stated otherwise, worms were harvested from Welsh ponies maintained in INRAE facilities. Infected ponies were given pyrantel pamoate (6.6 mg/kg liveweight, Strongid®, Zoetis, Malakoff, France) per os and faecal matter was harvested between 18 and 24 h after treatment to collect individual male worms. Upon collection, worms were bathed in PBS 1x to remove faecal debris, placed in microcentrifuge tubes left on ice for no more than 15 min before being put at -80C. For Hi-C data production, 70 worms were used and collected following the same procedure. In that case however, worms had their head cut for genotyping while their body was flash frozen.
Transcriptomic data were generated from whole worms. These worms were picked from faecal matter collected 18h after pyrantel treatment, briefly bathed in 1x PBS, placed in a microcentrifuge tube and flash frozen in liquid nitrogen.

To quantify the sex-ratio in cyathostomins released after pyrantel treatment, the faecal matter of 10 Welsh ponies was collected between 15 and 24 h after drug administration. On every occasion, the worms of each sex were counted from 200 g of faeces.

## DNA and RNA extraction protocols for genome assembly and annotation

DNA was extracted following a salting-out procedure to reduce the number of processing steps. Worms were gently crushed with a sterile pestle and incubated for 3 hours at 56°C in a lysis buffer adapted from a previous experiment on mosquitoes (Kingan et al., 2019). Following incubation, NaCl was added to the mix before precipitation was induced with isopropanol. After centrifugation (6250g at 4°C for 5 min), the DNA pellet was washed with 70% ethanol and resuspended in 30 µL of TE buffer. A second DNA library was produced using 2.5 µL native DNA from the same sample after whole genome amplification procedure following the REPLI-G® (Qiagen) manufacturer's recommendation.

Genome annotation was based on RNAseq data from male and female *C. nassatus*. To isolate *C. nassatus* from other cyathostomin species after worm collection, simultaneous DNA and RNA extraction (AllPrep® DNA/RNA mini kit, Qiagen) was performed on 94 worms. DNAs were used for Sanger sequencing of the ITS-2 and COI regions to subsequently isolate *C. nassatus* species. RNAs were stored at -80°C before three pools of four to five males and three pools of three to four females were made for sequencing.

## PacBio sequencing of a single *C. nassatus* worm

The REPLI-g amplified and non-amplified nativeDNA left from the same single male worm were subsequently processed at GeT-PlaGe core facility (INRAe Toulouse) to prepare two libraries according to the manufacturer's instructions "Procedure & Checklist - Preparing HiFi Libraries from Low DNA Input Using SMRTbell® Express Template Prep Kit 2.0". At each step, DNA was quantified using the Qubit dsDNA HS Assay Kit (Life Technologies). DNA purity was tested using a nanodrop (Thermofisher) and size distribution and degradation assessed using the Femto pulse Genomic DNA 165 kb Kit (Agilent). Purification steps were performed using AMPure PB beads (PacBio) and 1 µg of DNA was purified and then sheared at 13 kb using the Megaruptor1 system (Diagenode). The library was size-selected, using AMPure® PB beads (PacBio) to remove less than 3kb long templates. Using Binding kit 2.0 (primer V4, polymerase 2.0) and sequencing kit 2.0, the 12 kb library was sequenced onto 1 SMRTcell on Sequel2 instrument at 40 pM with a 2 hours pre-extension and a 30 hours movie.

A second library was prepared from the same DNA sample following the manufacturer's instructions "Procedure & Checklist - Preparing HiFi SMRTbell® Libraries from Ultra-Low DNA Input". At each step, DNA was quantified using the Qubit dsDNA HS Assay Kit (Life Technologies). DNA purity was tested using the nanodrop (Thermofisher) and size distribution and degradation assessed using the Femto pulse Genomic DNA 165 kb Kit (Agilent). Purification steps were performed using AMPure PB beads (PacBio). 13 ng of DNA was purified then sheared at 10 kb using the Megaruptor1 system (Diagenode). Using SMRTbell® gDNA Sample Amplification Kit, 5ng of DNA was amplified by 2 complementary PCR reactions (13 cycles). Then 500 ng of the library was size-selected, using a 6.5 kb cutoff on the BluePippin Size Selection system (Sage Science) with the "0.75% DF Marker S1 3-10 kb Improved Recovery" protocol. Using Binding kit 2.0 (primer V4, polymerase 2.0) and sequencing kit 2.0, the 9.5 kb library was sequenced onto one SMRTcell on Sequel2 instrument at 50 pM with a 2 hours pre-extension and a 30 hours movie. The low- and ultra-low input protocols yielded 6 and 26 Gbp respectively, corresponding to 615,427 and 2,780,682 CSS reads respectively and 10 and 44x depth of coverage.
## RNA sequencing of male and female *C. nassatus* worms

RNAseq was performed at the GeT-PlaGe core facility, INRAe Toulouse. RNA-seq libraries have been prepared according to Illumina's protocols using the Illumina TruSeq Stranded mRNA sample prep kit to analyse mRNA. Briefly, mRNA were selected using poly-T beads. Then, RNA were fragmented to generate double stranded cDNA and adaptors were ligated to be sequenced. 11 cycles of PCR were applied to amplify libraries. Library quality was assessed using a Fragment Analyser and libraries were quantified by QPCR using the Kapa Library Quantification Kit. RNA-seq experiments have been performed on an Illumina NovaSeq 6000 using a paired-end read length of 2x150 pb with the associated sequencing kits.

## Hi-C library preparation and sequencing

Hi-C library was constructed using the Arima-HiC kit (Arima, ref. A510008) and the Accel NGS 2S Plus DNA Library Kit (Swift Biosciences, ref. SW21024). Briefly, we pulverised complete worms in liquid nitrogen using a mortar and crosslinked the pulverised tissue using a 2% formaldehyde solution. After tissue lysis, we digested the crosslinked DNA according to the manufacturer's protocol. We repaired the digested DNA using biotinylated nucleotides and performed a ligation aiming for spatially proximal digested ends of DNA. We purified the proximally ligated DNA, sonicated it using a e220 focused-ultrasonicator (Covaris) and enriched the biotinylated fragments. Starting from enriched biotinylated fragments, we constructed a NGS library using the Accel-NGS 2S Plus DNA library kit (Swift Biosciences, SW21024) according to ARIMA's instruction, whereby fragments were first end-repaired before indexed adapters ligation. After purification, a small fraction of the indexed DNA was used to determine by qPCR the number of PCR cycles necessary for optimal amplification.

Based on this result, we performed 6 cycles of PCR amplification on the remaining indexed DNA. The size distribution of the resulting library was monitored using a Fragment Analyzer with the High Sensitivity NGS kit (Agilent Technologies, Santa Clara, CA, USA) and the library was quantified using microfluorimetry (Qubit dsDNA HS kit, Thermofischer scientific). The library was denatured with NaOH, neutralised with Tris-HCl, and diluted to 300 pM. Clustering and sequencing were performed on a NovaSeq 6000 (Illumina, San Diego, CA, USA) using the paired end 150 nt protocol on 1 lane of a flow cell SP. Image analyses and base calling were performed using the Miniseq Control Software and the Real-Time Analysis component (Illumina). Demultiplexing was performed using Illumina's conversion software (bcl2fastq 2.20.0.422). The quality of the raw data was assessed using FastQC v0.11.9 (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) from the Babraham Institute and the Illumina software SAV (Sequencing Analysis Viewer). FastqScreen (Wingett & Andrews, 2018) was used to identify potential contamination.

## Genome assembly, gene model prediction and annotation

HiFi data (30 Gb) were adapter-trimmed and filtered for low-quality reads. Genome assembly graphs were produced using HiFiasm (Cheng et al., 2021) v0.16.1 using default parameters. Assembly completeness was assessed using BUSCO (Seppey et al., 2019) v 5.2.2 (nematoda_odb10 gene set, n = 3,331) and metrics were estimated using Quast v-5.0.2 (Gurevich et al., 2013).

This first assembly was subsequently used as a back-bone for scaffolding with Hi-C data. These data were processed with the juicer pipeline (Durand et al., 2016a) v 1.5.7. The assembly was scaffolded with 3d-dna (Dudchenko et al., 2017) and manually corrected with juicebox (Durand et al., 2016b) (version 1.11.08).

The mitochondrial genome was assembled using the mitoHiFi v2.2 software (Uliano-Silva et al., 2023), feeding quality filtered HiFi reads and the *C. nassatus* reference mitogenome (Gao et al., 2017) as a backbone. The final mitogenome sequence was subsequently added to the assembly fasta file.
We also identified and removed the chimeric hifiasm assembled mitochondrion sequence (38 Kbp in length, matched to scaffold 57 with the minimap2 software (Li, 2018).

Repeat elements were identified using RepeatMasker (Tarailo-Graovac & Chen, 2009) v4.0.7, Dust (Morgulis et al., 2006) v1.0.0 and TRF (Benson, 1999) v4.09. To locate repeat regions, the *C. elegans* libraries and a *C. nassatus*-specific de novo repeat library built with RepeatModeler (Flynn et al., 2020) v1.0.11 were fed to RepeatMasker (Tarailo-Graovac & Chen, 2009). Bedtools v2.26.0 was used to aggregate repeated regions identified with the three tools and to soft mask the genome.

Gene models were built from RNAseq data generated from male and female samples. RNAseq reads were mapped onto the masked genome (repeatMasker v4.0.7) using the HiSat2 (Zhang et al., 2021) software to produce hints subsequently used for protein to genome alignment with exonerate. Ab initio gene structures were subsequently estimated with the BRAKER (Hoff et al., 2019) pipeline relying on both RNAseq data and available proteomes for other clade V species (*Ancylostoma ceylanicum*, *Haemonchus contortus*, *Caenorhabditis elegans*, *C. inopinata*, *C. remanei*, *C. tropicalis*) and more phylogenetically distant species (*Brugia malayi* and *Onchocerca volvulus* for clade III, and *Trichuris muris* for clade I parasitic nematodes). These gene predictions were fed to MAKER (Campbell et al., 2014) for final gene model prediction. Gene annotation was achieved with the EnTAP software (Hart et al., 2020) to assign homology information from the EggNOG (Hernández-Plaza et al., 2023) database and protein domain from PFAM (Mistry et al., 2021) and InterPro (Paysan-Lafosse et al., 2023) databases.

## Comparative genomic analyses

The synteny between the *C. nassatus* and *H. contortus* genomes was inferred after aligning both references against one another, using the PROMER software (Delcher et al., 2002) with the *-mum* option to recover the only exact matches being unique between the query and reference sequences. Hits were subsequently filtered with the delta-filter tool, setting the minimal length of a match at 500 amino acids. A custom perl script (Cotton et al., 2016) was subsequently used to draw the circus plot, showing links with 70% similarity and a minimal length of 10 Kbp. Ortholog groups were identified using the Orthofinder software v2.5.4 (Emms & Kelly, 2019), considering the subset of the most complete nematode genomes (N50 above 10 Mbp) present in WormBaseParasite v.17 and annotated with RNAseq data, including clade V parasitic (*H. contortus*) and free-living species (*C. briggsae*, *C. elegans*, *C. inopinata*, *C. nigoni*, *C. remanei*, *C. tropicalis* and *Pristionchus pacificus*), clade IV plant (*Bursaphelenchus xylophylus*, *Heterodera glycines*) and animal parasites (*Strongyloides ratti*), filarial species (*Brugia malayi*, *Onchocerca volvulus*) and *Trichuris muris* as a clade I representative. Because of the close phylogenetic relationship with ancylostomatidae, the *Ancylostoma ceylanicum* proteome was also considered. The species tree was built using the iTOL software (Letunic & Bork, 2021). Gene family evolution was performed with the CAFE v3 software (Han et al., 2013).

## Differential transcriptomic response of male and female *C. nassatus* after exposure to pyrantel

RNA-seq reads were mapped with the Salmon software v1.4 (Patro et al., 2017) with correction for GC-content, sequence-specific and position biases and the validateMappings option. Pseudo-counts were imported with the tximport package v1.18.0 and differential gene expression was estimated with the DESeq2 v1.30.1 package (Love et al., 2014). Out of the 22,568 expressed genes, 12,569 with at least ten counts in five replicates were considered for analysis. Female worm gene counts displayed a bimodal distribution. The cut-off splitting the two underlying distributions was determined with the peakPick package v0.11. Differential gene expression was applied to both gene populations, i.e. those with a median count below or above the identified cut-off. In any case, genes showing an absolute fold change above 2 and an adjusted P-value below 1% were deemed significant.

Gene Ontology enrichment was run using the topGO package v2.42 (Alexa et al., 2006) and processed with the GeneTonic package (Marini et al., 2021). GO terms with P-value below 1% were
considered as significant. To gain additional functional insights, one-to-one orthologs of the differentially expressed genes in *C. elegans* were considered for tissue or phenotype enrichment using the wormbase enrichment tool (considering a q-value cut-off of 1%). To investigate differences between male and female worms that would not be associated with the presence of egg in female utero, the subset of DE genes with one-to-one ortholog in *C. elegans* were considered to remove any gene known to be expressed in oocyte, the germ line or the embryo (data retrieved using the wormbase.org "simple mine" online tool). Genes with unknown expression patterns in *C. elegans* were not considered. While this subset remains artificial in nature, it offers a gene population whose differential expression between males and females was deemed independent from the presence of eggs in female worms.

## Spatio-temporal sampling of *C. nassatus* diversity using a pool-sequencing framework

The sampling scheme pursued two objectives driven by the understanding of the evolution of drug resistance in cyathostomin populations. First, a spatio-temporal sampling was applied to establish the degree of connectivity among populations of Western Europe and America, including isolates from western Ukraine (n = 1), Germany (Hannover region, n = 1), France (Normandy, n = 4; Nouzilly, n = 1) and American (Kentucky, n = 1). French samples consisted of the isolate used for genome assembly and four other populations were collected in Normandy (Calvados département, within a 30 Km radius; note that latitude and longitude have been voluntarily modified to ensure anonymity of the sampled stud farms). In that case, two pyrantel-resistant isolates were geographically matched with pyrantel susceptible isolates. Pyrantel-resistant worms were exposed to pyrantel embonate to first confirm their resistance potential (Strongid®, 6.6 mg/kg, Zoetis, France) 14 days after treatment, and finally recovered 24 hours after ivermectin treatment (Eqvalan pâte, 200 µg/kg, Boehringer Ingelheim animal health, France). Pyrantel sensitivity (measured as a Faecal Egg Count Reduction coefficient) was estimated using the eggCounts package (Wang et al., 2018) for the Normandy isolates, while previously published data were used for the reference (Boisseau et al., 2023) and American isolates (Scare et al., 2018).

For every population, single male worms were considered and genotyped for the ITS-2 and COI barcodes (Courtot et al., 2023) to ascertain their identity. DNA was extracted using the same salting-out protocol as before, quantified using a Qubit apparatus, and pooled together with an equi-molar contribution of each individual worm to the pool. The Ukrainian population consisted of 100 worms sampled after ivermectin treatment across Ukrainian operations and morphologically identified (Kuzmina, 2012). Their DNA was extracted en-masse and used for pool-sequencing.

A museum sample collected in 1899 by Prof. Arthur Looss in Egypt, fixed in 70% ethanol and preserved in the Swedish Museum of Natural History (Kuzmina & Holovachov, 2023) was also included in the study to compare pre-anthelmintic diversity to that observed in contemporary samples. In that case, worms were identified morphologically. To deal with this old material, the previously described extraction protocol (Gamba et al., 2016) was adapted as follows. Briefly, worms were digested in a lysis buffer for 3 hours at 56°C. Worm lysate was subsequently processed using the DNA clean-up kit (Macherey-Nagel) using the washing solution in a 4-to-1 excess to retain as much material as possible. Elution buffer was heated at 70°C for DNA recovery in a final volume of 20µL.

Due to the small size of the starting DNA material, the pool-seq library was constructed using the Accel NGS 2S Plus DNA Library Kit (Swift Biosciences, ref. SW21024) without any fragmentation step. Briefly, we end-repaired the fragments before indexed adapter ligation. We performed a 9 cycles PCR amplification on the indexed DNA. The size distribution of the resulting library was monitored using a Fragment Analyzer with the High Sensitivity NGS kit (Agilent Technologies, Santa Clara, CA, USA) and the library was quantified using the KAPA Library quantification kit (Roche, Basel, Switzerland).

The library was denatured with NaOH, neutralised with Tris-HCl, and diluted to 200 pM. Clustering and sequencing were performed on a NovaSeq 6000 (Illumina, San Diego, CA, USA) using the single
read 100 nt protocol on 1 lane of a flow cell SP.Image analyses and base calling were performed using the Miniseq Control Software and the Real-Time Analysis component (Illumina). Demultiplexing was performed using Illumina's conversion software (bcl2fastq 2.20.0.422). The quality of the raw data was assessed using FastQC (v0.11.9) from the Babraham Institute and the Illumina software SAV (Sequencing Analysis Viewer). FastqScreen was used to identify potential contamination.

## Pool-sequencing data processing and SNP calling

After adapter trimming using cutadapt v.3.4, genomic data were mapped onto the reference using bwa-mem2 v2.2, retaining properly mapped reads with quality phred score above 20 with samtools v1.9. Duplicate reads were identified and filtered using Picard v.2.20.7 (https://broadinstitute.github.io/picard/) and indel realignment was applied using GATK v3.8 (McKenna et al., 2010). The pre-processed data were subsequently treated with the Popoolation2 software to generate a sync file, considering a minimum quality phred score of 30 and removing regions falling within 5 bp from an indel.

Genome-wide nucleotide diversity was estimated within 5 kbp windows, with minimal read depth of 10 using npstats (Ferretti et al., 2013). The average nucleotide diversity estimates were used as an initial guess for SNP calling using snape-pooled v2.8 software (Raineri et al., 2012) applied on filtered mpileup files (-C 50 -q 20 -Q 30 -d 400). The sites with 90% support were retained as high-quality SNPs and the intersection between the snape-pooled SNPs and popoolation2 sites were further considered for analysis.

This approach was applied to the whole set of populations for the study of genome-wide diversity (n = 23,384,719 autosomal SNPs supported by one read across pool; n = 514,068, after filtering for a depth of coverage between 20 and 400x, a minor allele frequency of 5%) or restricted to the set of six modern populations with pyrantel efficacy data (n = 36,129,276 autosomal SNPs and 2,807,924 SNPs on the X chromosome) for the association study.

## Investigation of genome-wide diversity and pyrantel resistance architecture

Nucleotide diversity and Tajima's D statistics were estimated for 100-Kbp sliding windows with npstats v2.8 from the high quality SNP set provided as a bed file. Investigation of genome-wide diversity was subsequently performed on the set of SNPs shared across old and modern populations after filtering for sites with minimum coverage below 15, maximum coverage above 400 and minor allele frequency below 5% with the poolstat package v2.0.0 (Gautier et al., 2022). Private variants were selected from the set of unfiltered autosomal SNPs, isolating markers with a minimal depth of coverage of 10 in every pool but absent in all but the focal isolate. Reference allele frequencies were binned by 10% MAF for plotting. Population connectivity was investigated using a principal component analysis applied to folded allele frequencies using the pcadapt (Luu et al., 2017) R package v4.3.3 and the set of filtered SNPs (n = 1,346,424). Genetic differentiation (FST) and admixture statistics (f3 and f4) were estimated using the "anova" method of the poolstat v2.0.0 (Gautier et al., 2022).

To isolate genomic regions associated with pyrantel susceptibility, we regressed allele frequencies upon scaled population pyrantel efficacy while accounting for population genomic relatedness using the Baypass software v2.2 (Gautier, 2015). The modern population SNP set was filtered using the poolstat (Gautier et al., 2022) R package to retain sites with a pool coverage between 10 and 400, a minimal MAF of 5% and 4 reads to support an allele leaving 8,365,939 and 1,049,048 SNPs on the autosomes and X chromosome respectively. The analysis was run on 15 pseudo-replicates of 557,727 to 557,732 autosomal SNPs (69,936 to 69,936 for the X chromosome) and averaged from six independent runs (-npilot 15 -pilotlength 500 -burnin 2500). For the X chromosome, the analysis was run on the five isolates composed of males. Any variant with a Bayes factor above 35 was deemed decisive, while significance was decided upon a Bayes factor of 20. QTL regions were defined as the extreme edges of contiguous 1-Mbp windows with more than 10 decisive SNPs. Within these windows,
three strands of evidence were considered to call candidate genes. First, the known candidate genes for pyrantel resistance (unc-29, unc-38 and unc-63) were inspected for the presence of significant SNPs within their locus or in their vicinity. A second prior-informed approach consisted in inspecting homologs of *C. elegans* genes affected by aldicarb (https://wormbase.org/resources/molecule/WBMol:00003650#032--10-3, last accessed May 24th 2023) or levamisole (https://wormbase.org/resources/molecule/WBMol:00004019#032--10-3, last accessed May 24th 2023), two molecules affecting the effect of acetylcholine at the neuro-muscular junctions. Third, enrichment analysis was run using the *C. elegans* homologs with one-to-one or many-to-one orthology (e.g. four copies of *unc-29*) and the wormbase online enrichment tool (Lee et al., 2018). Fourth, the decisive SNPs within gene loci showing a correlation between allele frequency and the measured FECRT above 90% were regarded as strong candidates for field use, as their association would not be affected by population structure.

### C. elegans pyrantel-resistant line selection and RNAseq analysis

To generate pyrantel-resistant lines of *C. elegans*, a large population of wild-type N2 worms was produced from four 4th larval stage (L4) on 10 different NGM (Nematode Growth Medium) plates. All worms from this founder population were collected with M9 medium and pooled, before being transferred to new plates containing 20 µM pyrantel (IC₃₀ allowing 80% of worms reaching adulthood) or not, with 6 different replicates for each condition. Once the worms survived and reproduced (generally at least two generations per concentration) at a given concentration, they were transferred to a plate with higher pyrantel concentration with 20-µM increase until survival at 100 µM was reached. This resistance level was reached after 12 generations. Total RNA extractions were performed using Trizol reagent (Invitrogen, Carlsbad, CA, USA) and total RNA was isolated according to the manufacturer's recommendations. Briefly, worms from each condition were recovered with M9 and washed with M9, 3 times, before being resuspended in 750µL Trizol. They were then ground with the Precellys 24 Homogenizer (Bertin, France) using 10-15 glass beads (0.5 mm) and the following program: 6400 rpm, 10 s two times, twice for a total time of 40 s. After grinding, 250µL of Trizol were added to reach a final volume of 1 mL. A DNase treatment step was performed using rDNase from the NucleoSpin RNA XS kit (Macherey-Nagel, Germany) and the RNA concentrations were measured using a nanodrop spectrophotometer.

RNA library preparations, and sequencing reactions for both the projects were conducted at GENEWIZ Germany GmbH (Leipzig, Germany) as follows. RNA samples were quantified using Qubit 4.0 Fluorometer (Life Technologies, Carlsbad, CA, USA) and RNA integrity was checked with RNA Kit on Agilent 5600 Fragment Analyzer (Agilent Technologies, Palo Alto, CA, USA).

The strand-specific RNA sequencing library was prepared by using NEBNext Ultra II Directional RNA Library Prep Kit for Illumina following manufacturer's instructions (NEB, Ipswich, MA, USA). Briefly, the enriched RNAs were fragmented for 8 minutes at 94°C. First strand and second strand cDNA were subsequently synthesised. The second strand of cDNA was marked by incorporating dUTP during the synthesis. cDNA fragments were adenylated at 3'ends, and indexed adapters were ligated to cDNA fragments. Limited cycle PCR was used for library enrichment. The dUTP incorporated into the cDNA of the second strand enabled its specific degradation to maintain strand specificity. Sequencing libraries were validated using DNA Kit on the Agilent 5600 Fragment Analyzer (Agilent Technologies, Palo Alto, CA, USA), and quantified by using Qubit 4.0 Fluorometer (Invitrogen, Carlsbad, CA, USA).

The sequencing libraries were multiplexed and clustered onto a flow-cell on the Illumina NovaSeq instrument according to manufacturer's instructions. The samples were sequenced using a 2x150bp Paired End (PE) configuration. Image analysis and base calling were conducted by the NovaSeq Control Software (NCS). Raw sequence data (.bcl files) generated from Illumina NovaSeq was
converted into fastq files and de-multiplexed using Illumina bcl2fastq 2.20 software. One mis-match was allowed for index sequence identification.

Data were analysed using a similar framework as that applied to the sex-specific differential transcriptomic analysis using pseudo-mapping and performing differential analysis with the DESeq2 package (Love et al., 2014).

# Results

## The chromosome-shape assembly of *C. nassatus* defines an XO/XX karyotype

A reference genome was built from a single *C. nassatus* worm using PacBio long reads to build contigs which were scaffolded with Hi-C reads (Figure 1a, S1, Table S1) generated from a pool of 70 worms from the same isolate. This was part of a larger initiative that also produced long read data for five other species including *Coronocyclus labiatus*, *Cyathostomum catinatum*, *Cylicostephanus goldi*, *Cylicostephanus longibursatus* and *Cylicocyclus insigne* (see the Data subsection).

The first round of assembly using HiFi reads yielded a 666.9 Mbp sequence fragmented in 2,232 contigs after haplotype purging. Hi-C sequencing yielded 525 M read pairs, 55.41% of which being alignable to the assembly (Figure S1). Among these, 135,180,104 chromatin contacts were identified with an equal representation of inter- (n = 69,144,270) and intra-chromosomal (n = 66,035,834) levels and 36% of these contacts being found at 20 Kbp. Following Hi-C data curation, a chromosome-scale assembly of 514.7 Mbp was built (scaffold N50 = 91,661,033 bp, scaffold L50 = 3; contig N50 = 674,848 bp; Table 1, Figure 1a, Figure S1) with high degree of completeness as supported by the 85.8% complete BUSCOs found that slightly outperformed the *H. contortus* genome statistics of 82.7%. The difference between the two genome assemblies mostly owed to higher missing BUSCOs in the *H. contortus* assembly (9.6%, n = 300) relative to *C. nassatus* (6.4%, n = 203).

After RNA-seq informed gene model prediction, the genome was 94.1% gene complete (n = 2,841 single-copy and 188 duplicated BUSCOs) and 4.4% missing identifiers (n = 136). The annotation comprised 22,718 coding genes, with an average of 4,010.4 genes (range between 4,152 and 3,831 genes) sitting on the five larger chromosomes. Gene length (1,286 bp on average, range between 153 and 79,776 bp) was well correlated with the number of exon (Pearson's *r* = 0.76, P < 10⁻⁴) that could reach as high as 252 exons (*CNASG00031410*, chromosome I). Total interspersed repeats covered 50.07% of the assembly sequence (Figure 1, Table S2, Figure S2). Most of them (n = 656,357) were unclassified or otherwise falling in the LINEs (5.06% of sequence length), LTRs (1.89% of sequence length) or the Tc1/Mariner DNA families (n = 2,998 elements) (Figure S2). This pattern matches the previous description of the trichostrongylid *H. contortus* genome although fewer elements (n = 292,099) were unclassified in that case (Doyle et al., 2020).

Six major linkage groups defining nuclear chromosomes with evidence of telomeric repeats accounted for 97% of the assembly size. A 13,878 bp mitochondrial genome was resolved from single HiFi reads. Nuclear chromosome lengths ranged between 107.7 and 83.5 Mbp for five of them and a smaller sequence of 34.4 Mbp (scaffold 6). Using several strands of evidence, this chromosome likely defines the X chromosome whose dosage is associated with sex in other clade V nematode species. First, a 27.8x ± 0.9 drop in depth of coverage was observed for that chromosome in re-sequenced male populations (t₁ = -30.85, P < 10⁻⁴, Figure 1b), and its coverage was halved in pools of males relative to the pools containing males and females (βMales & females x x = 34.4 ± 1.8, t₁ = 19.04, P < 10⁻⁴, Figure 1b). Second, windowed estimates of nucleotide diversity in male worms were also significantly reduced on this chromosome relative to that observed on the rest of their genome (πchr1 to 5 = 0.008687 ± 3.070 x 10⁻⁶, πchr6 = 0.007662 ± 1.154 x 10⁻⁵, t₈₁₇₃₂ = 98.9, P < 10⁻⁴, Figure 1c). Third, this scaffold showed marked synteny with *H. contortus* X chromosome (42.5% of the 762 alignments involving scaffold 6; Figure S3). Last, the genes found on this chromosome were more often up-regulated in female worms relative to their male counterparts (54.4% of DE genes), whereas the opposite was true
on the other five chromosomes (35 to 49% of the DE genes, Figure 1a). This likely X chromosome showed a denser gene content (55.2 genes/Mbp vs. 43.08 genes/Mbp for the five larger chromosomes) and longer gene length (282.27 bp ± 37.68 more than on the other five chromosomes, t = 7.49, P < 10⁻⁴). Repeat elements also appeared relatively contained to chromosome arms (50% of repeat elements within the first and last 10 Mbp).

The estimated species tree (from 2,387 orthogroups shared across the selected species; n = 30,491 orthogroups in total), supported the intermediate phylogenetic position of *C. nassatus* between *Ancylostomatoidae* and *Trichostrongyloidae* (Figure 1d) as previously found for *C. goldi* (International Helminth Genomes, 2019).

In the lack of karyotyping data for any cyathostomin species, this fully resolved assembly defines a six-chromosome genome and a XX (female)/ XO (males) sex determination as seen in the trichostrongylid *H. contortus* or free-living caenorhabditids.

<figure><img src="image_4.png" /><figcaption>**Figure 1** - Structure of the *Cylicocyclus nassatus* genome</figcaption></figure>

a. The circos plot consists of six tracks for each chromosome assembled, with their averaged GC-content presented on the outer ring. The heatmap corresponds to the content in identified repetitive regions, with overall content, DNA, LINE, LTR and other families presented towards the center of the plot. The third layer represents the differentially expressed genes (represented as log fold change; dot size is proportional to this value) found between males and females taken as the reference level (positive values indicate over-expression in males). The innermost track shows the average depth of coverage for nine re-sequenced populations. b. Distribution of windowed depth of coverage across the six major scaffolds in male-only or female and male populations (bottom). c. Nucleotide diversity distribution for each chromosome. The orange and blue lines stand for mean π value of chromosome 6 or that of other chromosomes respectively. d. Phylogenetic tree (estimated from 2,387 protein sequences) showing gene family evolution (expanded or contracted; numbers in brackets) across parasitic nematode species of interest (coloured by clades).
Allele frequencies of known anthelmintic targets have changed over the last century

The genome-wide diversity of *C. nassatus* populations was investigated from nine populations sampled across five countries (Figure 2a), including an older population sampled from Egypt at the end of the 19th century. The windowed depth of coverage was 38.6x on average (Figure 2b) but lower in the older sample and in the Hi-C derived dataset from the reference population (18x in both cases). Despite DNA fragmentation of the old sample, the data yielded an even representation of the whole genome (Figure S6) and did not show evidence of strong degradation (Figure S7).

*Cylicocyclus nassatus* genetic diversity was high as illustrated by the 23,384,719 autosomal Single Nucleotide Polymorphisms identified or their average genome-wide nucleotide diversity estimates ranging from 6.06 x10⁻³ in the reference population to 9.07 x10⁻³ in the old isolate (Figure 2b, table S4). These estimates would be compatible with an average effective population size of 766,102 individuals.

<figure><img src="image_5.png" /><figcaption>**Figure 2** - Genome-wide diversity in modern and old *Cylicocyclus nassatus* isolates
a. Isolate repartition. b. Isolate sensitivity to pyrantel (black and red dashed lines correspond to original and targeted efficacies). Pool composition, windowed depth of coverage and genome-wide diversity coefficients distribution are represented for every isolate (OLD: 19th century Egyptian isolate; SBA, STB, STG and STL: four stud farms from Normandy, France; SDE: Hanover region, Germany; SKE: Kentucky, USA; STR: reference isolate, Nouzilly, France; UKR: Ukraine).</figcaption></figure>

The old Egyptian worms had both the highest nucleotide diversity estimates (0.087% difference with modern isolates, t = 37.19, P < 10⁻⁴) and the highest count of private variants, which would be compatible with a global loss of diversity in *C. nassatus* populations in the sampled isolates over the last century or higher diversity within African worms. Investigation of that latter hypothesis would require contemporary Egyptian worms that were not available for this study. The most extreme variance in allele frequencies between the 19th and 21st century isolates (n = 1,267 outlier windows of 100 SNPs; Figure S8) occurred on chromosomes 2 (n = 748, Figure S9) and 5 (n = 290, Figure S10).
Consistent outlier signals (n = 7 out of 8 comparison) were found over five genes coding for homologs of a D-aspartate oxidase (ddo-1 and ddo-2 in C. elegans) or lipase-domain containing genes (CNASG00085640, CNASG00085630). On chromosome 2, the signals were enriched over QTL regions associated with pyrantel resistance (see paragraph 5, Figure S9) and a single outlier window (between the old and Kentucky isolate) encompassed the β-tubulin isotype-1 locus (CNASG00077300) known to confer resistance to benzimidazoles (Kwa et al., 1993) (Figure S9). No signal was found over the β-tubulin isotype-2 coding gene on chromosome 5 (Figure S10). Along with these known drug resistance candidate genes, outlier differentiation occurred over homologs of nAChR subunit coding genes sitting on chromosome 2, namely *unc-63* and three homologs of *unc-29* (see paragraph 5, Figure S9).

## The connectivity of *C. nassatus* populations is compatible with ancient and pervasive admixture from Eastern Europe

Most SNPs were private to each isolate (between 51,580 and 375,369 SNPs; Figure 3a), with only 1,346,424 common SNPs matching the filtering criteria. Private variants segregated (98% of 1,801,440 SNPs) at low frequencies (<30%) especially in the Ukrainian isolate suggestive of important gene flow in this region (Figure 3a). On the contrary, the American isolate displayed 113 sites with frequency of 60% and higher (Figure 3a), reflecting accumulation of private variants permitted by their lower genetic connectivity to Europe, owing to both geographical distance and the isolation of their hosts since 1979.

The connectivity inferred from private variants was mirrored in the PCA on folded SFS, whereby the east-west gradient broadly defined the first component with the old Egyptian and contemporary American isolates positioned at each end (Figure 3b, S11). The second PCA axis was driven by the disconnection between these two samples from the modern European cluster (Figure 3b) and could reflect a temporal gradient, either owing to old variants present in 1899 in the Egyptian worms, or to the putative sampling of older ancestry (native or imported) in American worms. The close relationship entertained by modern European worms was evident on other components and would also suggest admixture between these isolates (Figure 3b, S11). Similarly, the genetic differentiation was also driven by the east-west gradient in a set of otherwise lowly differentiated isolates (global FST = 0.041, 95% c.i. = 0.04 - 0.042; Figure 3c). Last, while the negative Tajima's D found in every population would also suggest on-going admixture (Table S4), higher values were observed for the populations from closed herds (SKE, Tajima's D = -0.419 ± 0.69 for the Kentucky isolate, and -0.482 ± 0.86 for the reference STR isolate; Figure 2b). Despite these strands of evidence in favour of pervasive admixture, f3 coefficients - that measures the degree of resemblance between one population and two other source populations - were non-significant (|Z-score| > 12). This is compatible with either too little or too high admixture between populations, or genetic drift in the focal isolate erasing the signal since admixture.

In line with the latter hypothesis, we found significant evidence of shared ancestry between four French populations at lower allele frequencies (MAF<20%, n = 254,604 SNPs). In each case, the f3 with highest support consistently involved the Ukrainian isolate as a source proxy, in combination with other isolates of Normandy or German origins (Figure S12a). This would be compatible with a contribution from an Eastern-related population in relatively recent times that is absent or has been lost at higher frequencies. Of note, Ukraine retained most of the older variation as illustrated by the low genetic differentiation they entertained (FST = 0.022, Figure 3c) and the highest share of old variants (Figure 3d, n = 70,283 SNPs). This Eastern contribution had also the best support for admixture events detected in the old worms when considering high frequency variants (MAF >70% in every population, n = 351,262; Figure S12b). These strands of evidence would hence suggest that the eastern European contribution to *C. nassatus* ancestry is pervasive and precedes the 19th century. Consistent with this pattern, the most likely admixture graph (Figure 3e, S13) was compatible with a founder population closely related to the Ukrainian worm populations and from which all other isolates
derived. In that scenario, successive admixture events between past populations sharing ancestry
with the German and Kentucky isolates would have defined the European, American and Egyptian
isolates (Figure 3e).

<figure><img src="image_6.png" /><figcaption>**Figure 3** - Pervasive gene flow defines old and modern *C. nassatus* isolates
a. The number of private SNPs for a given isolate with allele frequency below or equal to the allele frequency displayed on the x-axis. b. Projection of the isolates on the first two components of a PCA on allele frequencies (n = 1,223,750 SNPs). c. Heatmap of genetic differentiation coefficients estimated from 1,223,750 SNPs spanning chromosomes 1 to 5 (the darker, the most differentiated). d. Heatmap of pairwise counts of SNPs observed only twice in the set of common SNPs. e. Admixture graph with best statistical support connecting the seven isolates from every sampled continent, with branch length given in drift units. Isolate abbreviations match figure 2 caption.</figcaption></figure>

Transcriptomic differences across sexes upon pyrantel exposure

Cyathostomin (reference isolate, STR) collection upon pyrantel treatment of their host indicated a disbalanced sex-ratio whereby three times more females than males were collected over the considered 9-hour window ($\beta_{sex} = 1.08 \pm 0.28$, $t_{66} = 3.917$, $P = 2 \times 10^{-4}$; Figure 4a, S14). This 20% male to female ratio was lower than previous reports ranging from 0.4 to 0.6 (Silva et al., 1999; da Silva Anjos & de Lurdes A Rodrigues, 2006; Sallé et al., 2018). Transcriptomic differences between males and females (three pools of four to five males and three pools of three to four females) collected after pyrantel treatment were quantified. Inspection of gene count distribution revealed a bimodal distribution in females (centred at a median count cut-off of 608, Figure 4b), applying over every chromosome (Figure S15) and reflecting the transcriptomic contrast between the female body and their reproductive tract and the developing eggs it harbours (Figure S15, tables S5 and S6). This was also mirrored in the 2,600 genes up-regulated in females (out of 5,860 differentially expressed between sexes, table S7) that defined significant enrichment for 120 biological processes, with an overwhelming contribution of nucleic acid metabolism, or chromatin organisation and segregation (Table S8)
compatible with egg division and production. In this respect, 138 genes were related to embryo development ending in birth or egg hatching (GO:0009792, *P* = 7.5 x 10⁻⁴, Figure 4c, Table S8). In sharp contrast, male worms expressed higher levels of transcripts related to ion homeostasis and pH regulation (n = 7 terms, Table S8), muscle contraction and metabolism (n = 4 terms, Table S8), neuronal system (n = 2 terms, Table S8), or pyrimidine-containing compound catabolic process (GO:0072529, *P* = 3.98 x 10⁻²).

To delineate somatic transcriptional differences between males and females not owing to the presence of eggs *in utero*, we selected *C. nassatus* DE genes with one-to-one ortholog in *C. elegans* that were not known to be expressed in any of the oocyte, germ line or embryonic stage of that model species. This process retained 148 genes up-regulated in males whose ontology enrichment highlighted neuropeptide signal pathway (GO:0007218, *P* = 8.9e-05), including orthologs of genes coding for FMRF like peptide (*flp-1*, -5, -17 and -21; Figure 4d, Table S9). These genes were also associated with an over-representation of locomotion-associated phenotypes or acetylcholinesterase inhibitor response variant in *C. elegans* (Table S10, Figure 4d). These data highlight the molecular substrates of sex differences in worms exposed to pyrantel.

## A handful of Quantitative Trait Loci (QTL) underpin pyrantel resistance

Transcriptomic data identified a plastic component of the response to pyrantel *in vivo* while the genes on which pyrantel selection acts upon may differ. To define the genetic landscape of pyrantel resistance, a genome-wide association analysis between SNP allele frequency and pyrantel resistance was run on six populations from France and Kentucky, USA (Figure 5, S16, Table S11).

<table><caption>**Table 1** - Quantitative Trait Loci (QTL) associated with pyrantel sensitivity with proposed candidate genes</caption><thead><tr><th>Chr.</th><th>QTL boundaries (bp)</th><th>Number of SNPs (total on chr.)</th><th>Candidate genes</th></tr></thead><tbody><tr><td>1</td><td>87424141 - 88111850</td><td>14 (15)</td><td>dkf-2*, erg-28*</td></tr><tr><td>2</td><td>46547030 - 53471131</td><td>150 (167)</td><td>ddo-3**, dys-1***, unc-29.4", unc-29.1*", unc-29.2*", unc-29.3*", gsnl-1**, rla-0**, yop-1**</td></tr><tr><td>2</td><td>81892469 - 82887338</td><td>8 (167)</td><td>chaf-1", endu-1**, fer-1", nduf-7**, rpl-24.1**</td></tr><tr><td>4</td><td>21218710 - 21548186</td><td>27 (31)</td><td>pgp-9</td></tr></tbody></table>

Bold names were covered by at least one decisive Single Nucleotide Polymorphism (SNP), others had one significant SNP and were at most 20 Kbp away from a decisive SNP; *: affected by aldicarb in *C. elegans*; **: enriched in the anal depressor muscle; #: physical interactants in *C. elegans*

This analysis found 222 SNPs with decisive association (Figure 5, S16,17) on the five autosomes (43.7% and 83.3% falling within or 5-Kbp away from a gene locus) while 43 SNPs on chromosome X reached the significant level (BF > 20 dB; 48.8% falling within gene loci).
<figure><img src="image_7.png" /><figcaption>**Figure 4** - Differential expression between pyrantel-exposed males and females informs on differential response between the two sexes. a. Mean male and female cyathostomin (reference isolate, Nouzilly, France) counts (52.6% C. *nassatus* as inferred from larval metabarcoding (Boisseau et al., 2023) collected between 15 and 24 hours after pyrantel treatment are represented along with the respective regression lines and confidence intervals. b. The distribution of the median log-transformed gene counts is represented for each sex, showing a bimodal distribution occurring in female C. *nassatus*. c. Female up-regulated genes are dominated by genes involved in embryo development. The absolute fold-change (relative to males) of the expression of the genes defining this GO term enrichment is plotted against their adjusted P-values. Panels match chromosomes and genes are annotated with their C. *elegans* orthologs when available. d. Genes up-regulated in males are enriched in locomotion-associated phenotypes in C. *elegans*. The size of each dot corresponds to the number of phenotype enrichments found in C. *elegans* and yellow dots indicate association with acetylcholinesterase inhibitor response variants.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 5** - Four Quantitative Trait Loci (QTL) regions define variation in pyrantel resistance. Left panels (a,d,f) represent the statistical support of the association (measured as a Bayes Factor, BF) between SNP allele frequency with pyrantel resistance on three chromosomes of interest (the red dashed line stands for the significance cut-off). The analysis was focused on contemporary isolates with known pyrantel sensitivity, i.e. four isolates from Normandy, France (two drug-resistant: STB and STL, and two drug-susceptible: SBA and STG), the reference drug-susceptible isolate (STR) and the drug-resistant isolate from Kentucky, USA (SKE). For each chromosome, the association between the raw allele frequency of the most likely candidate gene is plotted against the pyrantel sensitivity for every considered isolate (panels b & c for chromosome 2, e & g for chromosomes 1 and 4 respectively). h. Differential gene expression profile (between male and female worms exposed to pyrantel) for the candidate genes laying within identified QTL regions (colours match the respective chromosomes as in panels a, d, f). i. Significance and expression fold change for the genes differentially expressed between control and pyrantel-selected *Caenorhabditis elegans* lines.</figcaption></figure>
The autosomal associations were enriched within four major QTL regions located on three autosomes (Tables 1, S11) while chromosome 5 harboured another region with lower statistical support (Figure S16, Table S11). Chromosome 2 - where homologs of three known candidates for pyrantel resistance (*unc-63*, *unc-29* and *unc-38*) were found (Figure 5a) - harboured 75% of the decisive SNPs that were enriched within a broad 6-Mbp region (46.5 to 53.5 Mbp, Figure 5a) or a narrower peak centred at 82 Mbp (Table 1, Figure 5a). These QTL regions also exhibited strong distortion of their allele frequencies between modern pyrantel-resistant isolates (and one pyrantel susceptible isolate) and 19th century old worms (Figure S8). On the 5' end of this chromosome, two significant associations were found over a homolog of *T15B7.1*, while the *unc-63* homolog was 104,372 Kbp down-stream to these positions. However, none of the 154 SNPs found within the *unc-63* locus defined a robust association with pyrantel resistance (average BF = -6.032, max = 13) and a single SNP 2,038 bp upstream (12,303,713 bp) reached the significance level (20.6 dB). The same pattern applied over the *unc-38* locus (303 SNP, average BF = -6, max = 11.67 dB) that was 3.5 Mbp downstream from the closest associated SNP (82,887,338 bp, dB = 37.0). On the contrary, the broader QTL region harboured four homologs of the *unc-29* gene (Figure 5a), one of which (CNASG00064360) contains a single decisive SNP (Figure 5b) and two others reaching significant association levels (five SNPs within CNASG00064260 and one SNP over CNASG00064350). The genes (n = 1,636) found in this QTL region were associated with a significant enrichment for locomotion (GO:0040011, P = 3.7e-05, n = 25 genes vs. 10 expected) and axon development (GO:0040011, P = 4.3e-04, n = 24 genes vs. 11 expected) ontologies (first and third most significant ontology terms).

Along with *unc-29* copies, 31 other genes exhibited a SNP with decisive association with pyrantel resistance over that region and five others laid between 80 and 83 Mbp (Table S11). In that latter region, an homolog of the aldicarb resistant *fer-1* gene (Krajacic et al., 2013) (Figure 5c) and an homolog of *chaf-1* homolog were also covered with a decisive SNP (Table 1, Figure S17). Within the *fer-1* gene, two decisive and non-synonymous mutations (frameshift variants) were identified (Figure 5c, S17). Of note, this region also encompassed an homolog of the β-tubulin isotype-1 coding gene but the pervasive benzimidazole resistance across considered isolates rules out association with pyrantel sensitivity. To further prioritise candidate genes of functional interest, tissue enrichment analysis with the *C. elegans* homologs of this set of 37 genes found an over-representation of the anal depressor muscle in *C. elegans* (WBbt:0004292, fold-change = 8, q-value = 5.5e-04), defined by eight genes (*gsnl-1*, *rpl-24.1*, *rla-0*, *ddo-3*, *unc-29* and *yop-1* for the 47-53 Mbp QTL and *endu-1*, *nduf-7* for the downstream QTL region; Table 1).

Because pyrantel is thought to disrupt cholinergic signalling at the neuro-muscular junction, we also inspected homologs of *C. elegans* genes (n = 106) known to be affected by aldicarb (an acetylcholinesterase inhibitor) and levamisole (an nAChR agonist), other than *unc-29* and *fer-1* already covered by decisive SNPs. This approach identified two homologs of genes whose mutation confers aldicarb hypersensitivity (Table 1), i.e. *dys-1* (aldicarb hypersensitive) located on chromosome 2 (18,840 bp away from the closest decisive SNP at 49,510,860 and covered with five significant SNPs) and *erg-28* found on chromosome 1 (5,118 bp away from a decisive SNP at 87,813,389, Figure 5c,d). On that chromosome, the narrow QTL region harboured seven genes, out of which an homolog of *dkf-2* also represents a likely candidate due to its ability to modulate sensitivity to aldicarb in some *C. elegans* genetic backgrounds (Sieburth et al., 2005). As a last candidate of interest, a *pgp-9* homolog - whose expression is modulated by pyrantel in the ascarid *Parascaris univalens* (Martin et al., 2020) - harboured 19 decisive SNPs that defined the QTL on chromosome 4 (Table 1, Figure 5e,f, S17). Among these decisive SNPs, four defined major non-synonymous changes including a premature stop codon in exon 27 (SNP 21,242,866 bp) thereby defining likely explanatory mutations. For these two QTL regions however, genetic differentiation between modern isolates with known pyrantel sensitivity and the old worm isolate were not departing strongly from the genome average (Figure S18, 19) and may be indicative of soft sweep in these regions.
Across the decisive SNPs, the major reference allele (susceptible isolate) was not preferentially associated with higher pyrantel efficacy (52 out of the 97 SNPs; Figure 5b,d,f, S19). In addition, allele frequency of the decisive variants falling within gene loci (n = 97) displayed absolute Pearson's *r* between 64.3 to 98.3% but missed the 5% significance cut-off for 50 of them (*P*-value between 0.052 and 0.16) underscoring the effect of population structure in some cases. Of note, 28 of these SNPs were already segregating in the old worm isolate at a frequency compatible with their surmised sensitivity (Figure S17). In the lack of preferential association of rare alleles with sensitivity, it is likely that the remaining SNPs have arisen over the last century (and that they were not missed in the old worm isolate).

At the transcriptomic level, three quarters of the genes harbouring one decisive SNP (n = 37 out of 97; Figure 5g, S20) were differentially expressed between males and females collected upon pyrantel treatment, up-regulation in males being the rule in most cases (n = 27 genes) as seen for the most likely candidate genes (Figure 5g). Most candidate genes did not show significant differences in nucleotide diversity (FDR > 5%). However, slight deviations from the chromosome average were found over the *gsnl-1* (*t*<sub>27</sub> = -2.19, nominal *P* = 0.038) and *rpl-24.1* (*t*<sub>14</sub> = -2.32, nominal *P* = 0.035) in the pyrantel resistant isolates (Figure S20).

Altogether, the SNPs found within a short-list of five genes (*unc-29*, *chaf-1*, *endu-1*, *fer-1* and *pgp-9*) define a set of robust genetic markers not affected by population structure.

**Transcriptomic profiling of *Caenorhabditis elegans* lines selected for pyrantel resistance supports the polygenic nature of this trait and gives support for highlighted candidate genes**

To cross-validate the association signals found across field *C. nassatus* isolates, we compared the transcriptomic profiles of divergent *C. elegans* lines exposed to pyrantel or not for 12 generations. This data supported differential expression (adj. P<0.05) for ten out of the 14 genes identified within QTL regions (Figure 5i). However, a subset of 596 genes with higher fold-changes (Table S12) encompassed genes that defined significant enrichment for movement-related phenotypes (paralysis, q-value = 2.1e-4, or movement variant, q-value = 2.4e-06) and hypersensitivity to levamisole or nicotine (Table S13). The genes associated with hypersensitivity to these molecules were related to cuticle formation suggesting that permeability to the drug may be an indirect evasion strategy for the worms.

**Discussion**

This work delivers the first resolution of any cyathostomin genomes and provides the first genomic landscape of pyrantel resistance in any parasitic nematodes. We also gained insights into the temporal evolution of genome-wide diversity in *C. nassatus* over the past 150 years.

Whole-genome based evolutionary analysis found that *C. nassatus* was a closer relative of *Ancylostomatidae* than the *Trichostrongyloidea* as already reported using the *C. goldi* proteome (International Helminth Genomes, 2019). Similar to *A. ceylanicum* (Merchant et al., 2022), the *C. nassatus* genome harboured active endogenous viral elements that had distinct transcriptional patterns across sexes. It remains unknown however whether sexual metabolism affects their expression or if they play an active role in defining sexual traits. The similarity with *Ancylostomatidae* did not apply to their genome size that was about twice as large for *C. nassatus* as for *A. ceylanicum* (Schwarz et al., 2015) (*Ancylostomatidae*, 313 Mbp) but matched the 614-Mbp genome of *Teladorsagia circumcincta* (Hassan et al., 2023) (*Trichostrongylidae*). This difference in genome size matched a significant expansion of transposon elements that has been encountered across other helminth genomes (International Helminth Genomes, 2019). Of note, transposons were contained to the edge of the X chromosome arms, which might suggest limited colonisation and a recent differentiation of that chromosome, although the abundance of transposons in sex chromosome is a
poor correlate of evolutionary times in a wide range of species (Matsubara et al., 2006; Chalopin et al., 2015).

Access to past worm material offered an opportunity to track the evolution of genetic diversity in the past century for *C. nassatus*. The pattern of rare variants sharing would favour significant connectivity between old Egyptian worms and modern Western populations. However, this pattern was more pronounced with the Ukrainian isolates that seems to have contributed to every other isolates diversity in more recent times. This pattern was also found in the reconstructed admixture graph, although the surmised important extent of gene flow between parasite populations due to horse movement and the relatively limited size of the population set may have obscured this analysis (Patterson et al., 2012; Gautier et al., 2022). The origin of Cyathostominae remains unknown to date but they encompass a wide range of hosts including wild equids (Lichtenfels et al., 2008; Kuzmina et al., 2009, 2013, 2020; Tombak et al., 2021) and elephants (Chel et al., 2020). Their parasitic mode of life would hence suggest shared tracks of history with their hosts. In this respect, additional sampling of *C. nassatus* isolates along the road of horse domestication routes could confirm the contribution of contemporary horse expansion from the lower Volga-Don region (Librado et al., 2021) to the building of current *C. nassatus* populations. Access to horse coprolites would be key to sampling more ancient times as illustrated for middle-age dated *Trichuris trichiura* (Doyle et al., 2022b). Additional genetic profiling of isolates from wild African equids could also support the importance of horse domestication and management in the population structure of this species.

Similar to observations made between ancient and modern *T. trichiura* (Doyle et al., 2022b), lower nucleotide diversity was found in contemporary isolates relative to old North African *C. nassatus* individuals. In that respect, additional sampling effort from contemporary North African worms would help resolve the temporal and geographical contributions to the observed contrast in nucleotide diversity. Of note, the most important genetic differentiation occurred over pyrantel-resistance associated regions supporting the contribution of modern anthelmintic treatments in reshaping global parasite diversity. In this respect, the mapping of pyrantel resistance Quantitative Trait Loci (QTL) highlighted a few distinct regions among which a functional hub in the middle of chromosome 2 dominated. In this region, *unc-29* was the only known functional candidate already associated with pyrantel resistance. Our data however widen the previous picture with other mutations found within genes affecting cholinergic signalling either through regulation of receptors distribution at the synapse (*dys-1*), the regulation of cholinergic signalling including *fer-1* (Krajacic et al., 2013), *dkf-2* (Sieburth et al., 2005), and *chaf-1* (Gottschalk et al., 2005), or an *erg-28* homolog that is required for the function of the potassium SLO-1 channels (Cheung et al., 2020), the targets of emodepside (Welz et al., 2011). As such, the mechanistic landscape becomes more complex, as it remains to determine if a single mutation in any of these genes is sufficient or necessary to affect pyrantel resistance and if so, what mutation occurred first and what compensatory mutations may be needed. The broad QTL region identified would suggest at least that these candidates belong to a genomic hub that is less prone to recombination. Further investigation abrogating the function of the identified candidate genes within different genetic knock-out backgrounds, as applied for the study of aldicarb resistance (Sieburth et al., 2005), would contribute to disentangling the relative importance of each gene to the phenotype of interest while reconstructing the matrix of epistatic interactions between the proposed candidates.

Among the list of candidates identified in the QTL regions, *unc-29* was the only known functional candidate already associated with pyrantel resistance while the role of *pgp-9* was corroborated in *Parascaris univalens* exposed to pyrantel (Martin et al., 2020). However, *unc-38* seems to be ruled out in *C. nassatus* and the contribution of *unc-63* remains uncertain as a few SNP were associated with pyrantel resistance in the locus vicinity. While the associated SNPs certainly define a set of markers to be tested for monitoring of field cyathostomin isolates, the high number of SNPs reaching significance and the large number of differentially expressed genes in selected *C. elegans* lines would be compatible with a polygenic architecture for pyrantel resistance.
Of note, the close vicinity between the β-tubulin locus - that is known to affect benzimidazole resistance in cyathostomin (Kwa et al., 1995; Hodgkinson et al., 2008) or other nematode species (Sallé et al., 2019) - and a pyrantel-associated QTL may result in linked selection over that region by any of these two drugs. Under the assumption of conserved synteny between *C. nassatus* and *A. ceylanicum*, this genomic proximity may underpin the previously described synergistic effect of pyrantel and the pro-benzimidazole drug febantel in the management of *A. ceylanicum* (Hopkins & Gyr, 1991).

This detailed resolution sheds light on a complex network of genes whose function in the regulation of cholinergic signalling or the regulation of other drug targets opens novel perspectives for the use of drug combination in the field.

## Acknowledgements

The authors are greatly indebted to the breeders who granted access to their facilities and horses for sample collection in Normandy. Participation of T. Kuzmina in this study was partially supported by the EU through the Recovery and Resilience Plan for Slovakia under project No.09I03-03-V01-00015.

Preprint version 2 of this article has been peer-reviewed and recommended by Peer Community In Genomics (https://doi.org/10.24072/pci.genomics.100290; Pollet, 2025).

## Data, scripts, code, and supplementary information availability

The raw sequencing data have been submitted to the European Nucleotide Archive under study accession PRJEB63274. The HiFi reads and Hi-C data generated for *C. nassatus* are deposited with accession ERS15765218 and ERS15765233, while HiFi data generated for *Coronocyclus labiatus*, *Cyathostomum catinatum*, *Cylicostephanus goldi*, *Cylicostephanus longibursatus* and *Cylicocyclus insigne* correspond to accessions ERS15970850, ERS15970852, ERS15970851, ERS15978829, ERS15970849 respectively.

Scripts, code, and supplementary information are available online (https://doi.org/10.57745/1FB3SD, The cyathostomin genomics consortium, 2023).

## Conflict of interest disclosure

The authors declare that they comply with the PCI rule of having no financial conflicts of interest in relation to the content of the article.

## Funding

This project was funded by the French horse and horse riding institute (IFCE) and the Fonds Éperon. MGX and Get-PlaGe (https://doi.org/10.15454/1.5572370921303193E12) acknowledge financial support from France Génomique National infrastructure, funded as part of "Investissement d'Avenir" program managed by Agence Nationale pour la Recherche (contract ANR-10-INBS-09).

## References

Alexa A, Rahnenführer J, Lengauer T (2006) Improved scoring of functional groups from gene expression data by decorrelating GO graph structure. *Bioinformatics*, **22**, 1600–1607.
https://doi.org/10.1093/bioinformatics/btl140

Antonopoulos A, Doyle SR, Bartley DJ, Morrison AA, Kaplan R, Howell S, Neveu C, Busin V, Devaney E, Laing R (2022) Allele specific PCR for a major marker of levamisole resistance in *Haemonchus*
*contortus*. *International Journal for Parasitology. Drugs and Drug Resistance*, **20**, 17–26.
https://doi.org/10.1016/j.ijpddr.2022.08.001

Bellaw JL, Nielsen MK (2020) Meta-analysis of cyathostomin species-specific prevalence and relative abundance in domestic horses from 1975-2020: emphasis on geographical region and specimen collection method. *Parasites & Vectors*, **13**, 509. https://doi.org/10.1186/s13071-020-04396-5

Benson G (1999) Tandem repeats finder: a program to analyze DNA sequences. *Nucleic Acids Research*, **27**, 573–580. https://doi.org/10.1093/nar/27.2.573

Blanchard A, Guegnard F, Charvet CL, Crisford A, Courtot E, Sauve C, Harmache A, Duguet T, O'Connor V, Castagnone-Sereno P, Reaves B, Wolstenholme AJ, Beech RN, Holden-Dye L, Neveu C (2018) Deciphering the molecular determinants of cholinergic anthelmintic sensitivity in nematodes: When novel functional validation approaches highlight major differences between the model *Caenorhabditis elegans* and parasitic species. *PLoS Pathog*, **14**, e1006996.
https://doi.org/10.1371/journal.ppat.1006996

Boisseau M, Dhorne-Pollet S, Bars-Cortina D, Courtot É, Serreau D, Annonay G, Lluch J, Gesbert A, Reigner F, Sallé G, Mach N (2023) Species interactions, stability, and resilience of the gut microbiota - helminth assemblage in horses. *iScience*, **106044**.
https://doi.org/10.1016/j.isci.2023.106044

Bradley M, Taylor R, Jacobson J, Guex M, Hopkins A, Jensen J, Leonard L, Waltz J, Kuykens L, Sow PS, Madeja U-D, Hida T, Ole-Moiyoi K, King J, Argaw D, Mohamed J, Polo MR, Yajima A, Ottesen E (2021) Medicine donation programmes supporting the global drive to end the burden of neglected tropical diseases. *Transactions of The Royal Society of Tropical Medicine and Hygiene*, **115**, 136–144. https://doi.org/10.1093/trstmh/traa167

Bucknell DG, Gasser RB, Beveridge I (1995) The prevalence and epidemiology of gastrointestinal parasites of horses in Victoria, Australia. *International Journal for Parasitology*, **25**, 711–724.
https://doi.org/10.1016/0020-7519(94)00214-9

Campbell MS, Holt C, Moore B, Yandell M (2014) Genome Annotation and Curation Using MAKER and MAKER-P. *Current Protocols in Bioinformatics*, **48**.
https://doi.org/10.1002/0471250953.bi0411s48

Chalopin D, Volff J-N, Galiana D, Anderson JL, Schartl M (2015) Transposable elements and early evolution of sex chromosomes in fish. *Chromosome Research*, **23**, 545–560.
https://doi.org/10.1007/s10577-015-9490-8

Chel HM, Iwaki T, Hmoon MM, Thaw YN, Soe NC, Win SY, Bawm S, Htun LL, Win MM, Oo ZM, Masum MA, Ichii O, Nakao R, Nonaka N, Katakura K (2020) Morphological and molecular identification of cyathostomine gastrointestinal nematodes of Murshidia and Quilonia species from Asian elephants in Myanmar. *International Journal for Parasitology. Parasites and Wildlife*, **11**, 294–301. https://doi.org/10.1016/j.ijppaw.2020.03.005

Cheng H, Concepcion GT, Feng X, Zhang H, Li H (2021) Haplotype-resolved de novo assembly using phased assembly graphs with hifiasm. *Nature Methods*, **18**, 170–175.
https://doi.org/10.1038/s41592-020-01056-5

Cheung TP, Choe J-Y, Richmond JE, Kim H (2020) BK channel density is regulated by endoplasmic reticulum associated degradation and influenced by the SKN-1A/NRF1 transcription factor (AC Hart, Ed.). *PLOS Genetics*, **16**, e1008829. https://doi.org/10.1371/journal.pgen.1008829

Collobert-Laugier C, Hoste H, Sevin C, Dorchies P (2002) Prevalence, abundance and site distribution of equine small strongyles in Normandy, France. *Veterinary Parasitology*, **110**, 77–83.
https://doi.org/10.1016/S0304-4017(02)00328-X

Cotton JA, Bennuru S, Grote A, Harsha B, Tracey A, Beech R, Doyle SR, Dunn M, Hotopp JC, Holroyd N, Kikuchi T, Lambert O, Mhashilkar A, Mutowo P, Nursimulu N, Ribeiro JM, Rogers MB, Stanley E, Swapna LS, Tsai IJ, Unnasch TR, Voronin D, Parkinson J, Nutman TB, Ghedin E, Berriman M,
Lustigman S (2016) The genome of *Onchocerca volvulus*, agent of river blindness. *Nat Microbiol*, **2**, 16216. https://doi.org/10.1038/nmicrobiol.2016.216

Courtot É, Boisseau M, Dhorne-Pollet S, Serreau D, Gesbert A, Reigner F, Basiaga M, Kuzmina T, Lluch J, Annonay G, Kuchly C, Diekmann I, Krücken J, Von Samson-Himmelstjerna G, Mach N, Sallé G (2023) Comparison of two molecular barcodes for the study of equine strongylid communities with amplicon sequencing. *PeerJ*, **11**, e15124. https://doi.org/10.7717/peerj.15124

Courtot E, Charvet CL, Beech RN, Harmache A, Wolstenholme AJ, Holden-Dye L, O'Connor V, Peineau N, Woods DJ, Neveu C (2015) Functional Characterization of a Novel Class of Morantel-Sensitive Acetylcholine Receptors in Nematodes (RV Aroian, Ed.). *PLOS Pathogens*, **11**, e1005267. https://doi.org/10.1371/journal.ppat.1005267

Cwiklinski K, Merga JY, Lake SL, Hartley C, Matthews JB, Paterson S, Hodgkinson JE (2013) Transcriptome analysis of a parasitic clade V nematode: comparative analysis of potential molecular anthelmintic targets in *Cylicostephanus goldi*. *Int J Parasitol*, **43**, 917-27. https://doi.org/10.1016/j.ijpara.2013.06.010

DALYs, Hale Collaborators (2017) Global, regional, and national disability-adjusted life-years (DALYs) for 333 diseases and injuries and healthy life expectancy (HALE) for 195 countries and territories, 1990-2016: a systematic analysis for the Global Burden of Disease Study 2016. *Lancet*, **390**, 1260-1344. https://doi.org/10.1016/S0140-6736(17)32130-X

Delcher AL, Phillippy A, Carlton J, Salzberg SL (2002) Fast algorithms for large-scale genome alignment and comparison. *Nucleic Acids Research*, **30**, 2478-2483. https://doi.org/10.1093/nar/30.11.2478

Doyle SR (2022) Improving helminth genome resources in the post-genomic era. *Trends in Parasitology*, **38**, 831-840. https://doi.org/10.1016/j.pt.2022.06.002

Doyle SR, Illingworth CJR, Laing R, Bartley DJ, Redman E, Martinelli A, Holroyd N, Morrison AA, Rezansoff A, Tracey A, Devaney E, Berriman M, Sargison N, Cotton JA, Gilleard JS (2019) Population genomic and evolutionary modelling analyses reveal a single major QTL for ivermectin drug resistance in the pathogenic nematode, *Haemonchus contortus*. *BMC Genomics*, **20**, 218. https://doi.org/10.1186/s12864-019-5592-6

Doyle SR, Laing R, Bartley D, Morrison A, Holroyd N, Maitland K, Antonopoulos A, Chaudhry U, Flis I, Howell S, McIntyre J, Gilleard JS, Tait A, Mable B, Kaplan R, Sargison N, Britton C, Berriman M, Devaney E, Cotton JA (2022a) Genomic landscape of drug response reveals mediators of anthelmintic resistance. *Cell Reports*, **41**, 111522. https://doi.org/10.1016/j.celrep.2022.111522

Doyle SR, Søe MJ, Nejsum P, Betson M, Cooper PJ, Peng L, Zhu X-Q, Sanchez A, Matamoros G, Sandoval GAF, Cutillas C, Tchuenté L-AT, Mekonnen Z, Ame SM, Namwanje H, Levecke B, Berriman M, Fredensborg BL, Kapel CMO (2022b) Population genomics of ancient and modern *Trichuris trichiura*. *Nature Communications*, **13**, 3888. https://doi.org/10.1038/s41467-022-31487-x

Doyle SR, Tracey A, Laing R, Holroyd N, Bartley D, Bazant W, Beasley H, Beech R, Britton C, Brooks K, Chaudhry U, Maitland K, Martinelli A, Noonan JD, Paulini M, Quail MA, Redman E, Rodgers FH, Sallé G, Shabbir MZ, Sankaranarayanan G, Wit J, Howe KL, Sargison N, Devaney E, Berriman M, Gilleard JS, Cotton JA (2020) Genomic and transcriptomic variation defines the chromosome-scale assembly of *Haemonchus contortus*, a model gastrointestinal worm. *Communications Biology*, **3**, 656. https://doi.org/10.1038/s42003-020-01377-3

Dudchenko O, Batra SS, Omer AD, Nyquist SK, Hoeger M, Durand NC, Shamim MS, Machol I, Lander ES, Aiden AP, Aiden EL (2017) De novo assembly of the *Aedes aegypti* genome using Hi-C yields chromosome-length scaffolds. *Science* (*New York*, *N.Y.*), **356**, 92-95. https://doi.org/10.1126/science.aal3327
Durand NC, Shamim MS, Machol I, Rao SSP, Huntley MH, Lander ES, Aiden EL (2016a) Juicer Provides a One-Click System for Analyzing Loop-Resolution Hi-C Experiments. *Cell Systems*, **3**, 95-98. https://doi.org/10.1016/j.cels.2016.07.002

Durand NC, Robinson JT, Shamim MS, Machol I, Mesirov JP, Lander ES, Aiden EL (2016b) Juicebox Provides a Visualization System for Hi-C Contact Maps with Unlimited Zoom. *Cell Systems*, **3**, 99-101. https://doi.org/10.1016/j.cels.2015.07.012

Durrant C, Thiele EA, Holroyd N, Doyle SR, Sallé G, Tracey A, Sankaranarayanan G, Lotkowska ME, Bennett HM, Huckvale T, Abdellah Z, Tchindebet O, Wossen M, Logora MSY, Coulibaly CO, Weiss A, Schulte-Hostedde AI, Foster JM, Cleveland CA, Yabsley MJ, Ruiz-Tiben E, Berriman M, Eberhard ML, Cotton JA (2020) Population genomic evidence that human and animal infections in Africa come from the same populations of *Dracunculus medinensis*. *PLOS Neglected Tropical Diseases*, **14**, e0008623. https://doi.org/10.1371/journal.pntd.0008623

Emms DM, Kelly S (2019) OrthoFinder: phylogenetic orthology inference for comparative genomics. *Genome Biology*, **20**, 238. https://doi.org/10.1186/s13059-019-1832-y

Ferretti L, Ramos-Onsins SE, Pérez-Enciso M (2013) Population genomics from pool sequencing. *Molecular Ecology*, **22**, 5561-5576. https://doi.org/10.1111/mec.12522

Flynn JM, Hubley R, Goubert C, Rosen J, Clark AG, Feschotte C, Smit AF (2020) RepeatModeler2 for automated genomic discovery of transposable element families. *Proceedings of the National Academy of Sciences*, **117**, 9451-9457. https://doi.org/10.1073/pnas.1921046117

Gamba C, Hanghøj K, Gaunitz C, Alfarhan AH, Alquraishi SA, Al-Rasheid KAS, Bradley DG, Orlando L (2016) Comparing the performance of three ancient DNA extraction methods for high-throughput sequencing. *Molecular Ecology Resources*, **16**, 459-469. https://doi.org/10.1111/1755-0998.12470

Gandasegui J, Onwuchekwa C, Krolewiecki AJ, Doyle SR, Pullan RL, Enbiele W, Kepha S, Hatherell HA, van Lieshout L, Cambra-Pellejà M, Escola V, Muñoz J (2022) Ivermectin and albendazole coadministration: opportunities for strongyloidiasis control. *The Lancet. Infectious Diseases*, **22**, e341-e347. https://doi.org/10.1016/S1473-3099(22)00369-3

Gao Y, Qiu J-H, Zhang B-B, Su X, Fu X, Yue D-M, Wang C-R (2017) Complete mitochondrial genome of parasitic nematode *Cylicocyclus nassatus* and comparative analyses with *Cylicocyclus insigne*. *Experimental Parasitology*, **172**, 18-22. https://doi.org/10.1016/j.exppara.2016.11.005

Gautier M (2015) Genome-Wide Scan for Adaptive Divergence and Association with Population-Specific Covariates. *Genetics*, **201**, 1555-79. https://doi.org/10.1534/genetics.115.181453

Gautier M, Vitalis R, Flori L, Estoup A (2022) *f*-Statistics estimation and admixture graph construction with Pool-Seq or allele count data using the R package *poolfstat*. *Molecular Ecology Resources*, **22**, 1394-1416. https://doi.org/10.1111/1755-0998.13557

Gottschalk A, Almedom RB, Schedletzky T, Anderson SD, Yates JR, Schafer WR (2005) Identification and characterization of novel nicotinic receptor-associated proteins in Caenorhabditis elegans. *The EMBO Journal*, **24**, 2566-2578. https://doi.org/10.1038/sj.emboj.7600741

Gurevich A, Saveliev V, Vyahhi N, Tesler G (2013) QUAST: quality assessment tool for genome assemblies. *Bioinformatics* (Oxford, England), **29**, 1072-1075. https://doi.org/10.1093/bioinformatics/btt086

Han MV, Thomas GWC, Lugo-Martinez J, Hahn MW (2013) Estimating gene gain and loss rates in the presence of error in genome assembly and annotation using CAFE 3. *Molecular Biology and Evolution*, **30**, 1987-1997. https://doi.org/10.1093/molbev/mst100

Hart AJ, Ginzburg S, Xu M (Sam), Fisher CR, Rahmatpour N, Mitton JB, Paul R, Wegrzyn JL (2020) EnTAP: Bringing faster and smarter functional annotation to non-model eukaryotic transcriptomes. *Molecular Ecology Resources*, **20**, 591-604. https://doi.org/10.1111/1755-0998.13106

Hassan SU, Chua EG, Paz EA, Tay CY, Greeff JC, Palmer DG, Dudchenko O, Aiden EL, Martin GB, Kaur P (2023) Chromosome-length genome assembly of *Teladorsagia circumcincta* - a globally
important helminth parasite in livestock. *BMC Genomics*, **24**, 74. https://doi.org/10.1186/s12864-023-09172-0

Herd RP (1990) The changing world of worms: the rise of the cyathostomes and the decline of *Strongylus vulgaris*. , **12**, 732–734.

Hernández-Plaza A, Szklarczyk D, Botas J, Cantalapiedra CP, Giner-Lamia J, Mende DR, Kirsch R, Rattei T, Letunic I, Jensen LJ, Bork P, von Mering C, Huerta-Cepas J (2023) eggNOG 6.0: enabling comparative genomics across 12 535 organisms. *Nucleic Acids Research*, **51**, D389–D394.
https://doi.org/10.1093/nar/gkac1022

Hodgkinson JE, Clark HJ, Kaplan RM, Lake SL, Matthews JB (2008) The role of polymorphisms at beta tubulin isotype 1 codons 167 and 200 in benzimidazole resistance in cyathostomins. *Int J Parasitol*, **38**, 1149–60. https://doi.org/10.1016/j.ijpara.2008.02.001

Hoff KJ, Lomsadze A, Borodovsky M, Stanke M (2019) Whole-Genome Annotation with BRAKER. In: *Gene Prediction* Methods in Molecular Biology. (ed Kollmar M), pp. 65–95. Springer New York, New York, NY. https://doi.org/10.1007/978-1-4939-9173-0_5

Hopkins TJ, Gyr P (1991) Synergism of a combination of febantel and pyrantel embonate against *Ancylostoma caninum* on dogs. *Vet Med Rev*, **61**, 3–9.

International Helminth Genomes C (2019) Comparative genomics of the major parasitic worms. *Nat Genet*, **51**, 163–174. https://doi.org/10.1038/s41588-018-0262-1

Jenkins E, Backwell A-L, Bellaw J, Colpitts J, Liboiron A, McRuer D, Medill S, Parker S, Shury T, Smith M, Tschritter C, Wagner B, Poissant J, McLoughlin P (2020) Not playing by the rules: Unusual patterns in the epidemiology of parasites in a natural population of feral horses (Equus caballus) on Sable Island, Canada. *International Journal for Parasitology: Parasites and Wildlife*, **11**, 183–190. https://doi.org/10.1016/j.ijppaw.2020.02.002

Jürgenschellert L, Krücken J, Bousquet E, Bartz J, Heyer N, Nielsen MK, von Samson-Himmelstjerna G (2022) Occurrence of Strongylid Nematode Parasites on Horse Farms in Berlin and Brandenburg, Germany, With High Seroprevalence of Strongylus vulgaris Infection. *Frontiers in Veterinary Science*, **9**, 892920. https://doi.org/10.3389/fvets.2022.892920

Kaplan RM, Vidyashankar AN (2012) An inconvenient truth: Global worming and anthelmintic resistance. *Vet Parasitol*, **186**, 70–78. https://doi.org/10.1016/j.vetpar.2011.11.048

Kingan SB, Heaton H, Cudini J, Lambert CC, Baybayan P, Galvin BD, Durbin R, Korlach J, Lawniczak MKN (2019) A High-Quality De novo Genome Assembly from a Single Mosquito Using PacBio Sequencing. *Genes*, **10**, 62. https://doi.org/10.3390/genes10010062

Kopp SR, Coleman GT, Traub RJ, McCarthy JS, Kotze AC (2009) Acetylcholine receptor subunit genes from *Ancylostoma caninum*: Altered transcription patterns associated with pyrantel resistance. *International Journal for Parasitology*, **39**, 435–441. https://doi.org/10.1016/j.ijpara.2008.08.005

Kopp SR, Kotze AC, McCarthy JS, Traub RJ, Coleman GT (2008) Pyrantel in small animal medicine: 30 years on. *The Veterinary Journal*, **178**, 177–184. https://doi.org/10.1016/j.tvjl.2007.06.021

Krajacic P, Pistilli EE, Tanis JE, Khurana TS, Lamitina ST (2013) FER-1/Dysferlin promotes cholinergic signaling at the neuromuscular junction in *C. elegans* and mice. *Biology Open*, **2**, 1245–1252.
https://doi.org/10.1242/bio.20135637

Kuzmina T (2012) Analysis of Regional Peculiarities of Strongylid (Nematoda, Strongylidae) Biodiversity in Domestic Horses in Ukraine. *Vestnik Zoologii*, **46**, e-7-e-15.

Kuzmina TA, Dzeverin I, Kharchenko VA (2016) Strongylids in domestic horses: Influence of horse age, breed and deworming programs on the strongyle parasite community. *Vet Parasitol*, **227**, 56–63. https://doi.org/10.1016/j.vetpar.2016.07.024

Kuzmina T, Holovachov O (2023) Equine Strongylidae and other parasitic nematodes described by Arthur Looss during 1895–1911 in the collections of the Swedish Museum of Natural History. *Zootaxa*, **5227**, 151–193. https://doi.org/10.11646/zootaxa.5227.2.1
Kuzmina T, Kharchenko V, Zvegintsova N, Zhang L, Liu J (2013) Strongylids (Nematoda: Strongylidae) in two zebra species from the “Askania-Nova” Biosphere Reserve, Ukraine: biodiversity and parasite community structure. *Helminthologia*, **50**, 172–180.
https://doi.org/10.2478/s11687-013-0128-0

Kuzmina T, Kuzmin YI (2008) The strongylid community of working donkeys (*Equus asinus* L.) in Ukraine. *Vestnik Zoologii*, **42**, 99–104.

Kuzmina TA, Zvegintsova NS, Yasynetska NI, Kharchenko VA (2020) Anthelmintic resistance in strongylids (Nematoda: Strongylidae) parasitizing wild and domestic equids in the Askania Nova Biosphere Reserve, Ukraine. *Annals of Parasitology*, **66**, 49–60.

Kuzmina T, Zvegintsova N, Zharkikh T (2009) Strongylid community structure of the Przewalski's horses (*Equus ferus przewalskii*) from the Biosphere Reserve "Askania-Nova", Ukraine. *Vestnik Zoologii*, **43**, 209–215.

Kwa MS, Kooyman FN, Boersema JH, Roos MH (1993) Effect of selection for benzimidazole resistance in *Haemonchus contortus* on beta-tubulin isotype 1 and isotype 2 genes. *Biochem Biophys Res Commun*, **191**, 413-9. https://doi.org/10.1006/bbrc.1993.1233

Kwa MS, Veenstra JG, Van Dijk M, Roos MH (1995) Beta-tubulin genes from the parasitic nematode *Haemonchus contortus* modulate drug resistance in *Caenorhabditis elegans*. *J Mol Biol*, **246**, 500–10.
https://doi.org/10.1006/jmbi.1994.0102

Laing R, Doyle SR, McIntyre J, Maitland K, Morrison A, Bartley DJ, Kaplan R, Chaudhry U, Sargison N, Tait A, Cotton JA, Britton C, Devaney E (2022) Transcriptomic analyses implicate neuronal plasticity and chloride homeostasis in ivermectin resistance and response to treatment in a parasitic nematode. *PLOS Pathogens*, **18**, e1010545. https://doi.org/10.1371/journal.ppat.1010545

Laing R, Gillan V, Devaney E (2017) Ivermectin - Old Drug, New Tricks? *Trends Parasitol*, **33**, 463–472.
https://doi.org/10.1016/j.pt.2017.02.004

Lee RYN, Howe KL, Harris TW, Arnaboldi V, Cain S, Chan J, Chen WJ, Davis P, Gao S, Grove C, Kishore R, Muller HM, Nakamura C, Nuin P, Paulini M, Raciti D, Rodgers F, Russell M, Schindelman G, Tuli MA, Van Auken K, Wang Q, Williams G, Wright A, Yook K, Berriman M, Kersey P, Schedl T, Stein L, Sternberg PW (2018) WormBase 2017: molting into a new stage. *Nucleic Acids Res*, **46**, D869–D874.
https://doi.org/10.1093/nar/gkx998

Letunic I, Bork P (2021) Interactive Tree Of Life (iTOL) v5: an online tool for phylogenetic tree display and annotation. *Nucleic Acids Research*, **49**, W293–W296.
https://doi.org/10.1093/nar/gkab301

Li H (2018) Minimap2: pairwise alignment for nucleotide sequences (I Birol, Ed.). *Bioinformatics*, **34**, 3094–3100.
https://doi.org/10.1093/bioinformatics/bty191

Librado P, Khan N, Fages A, Kusliy MA, Suchan T, Tonasso-Calvière L, Schiavinato S, Alioglu D, Fromentier A, Perdereau A, Aury J-M, Gaunitz C, Chauvey L, Seguin-Orlando A, Der Sarkissian C, Southon J, Shapiro B, Tishkin AA, Kovalev AA, Alquraishi S, Alfarhan AH, Al-Rasheid KAS, Seregély T, Klassen L, Iversen R, Bignon-Lau O, Bodu P, Olive M, Castel J-C, Boudadi-Maligne M, Alvarez N, Germonpré M, Moskal-del Hoyo M, Wilczyński J, Pospuła S, Lasota-Kuś A, Tunia K, Nowak M, Rannamäe E, Saarma U, Boeskorov G, Lõugas L, Kyselý R, Peške L, Bălășescu A, Dumitrașcu V, Dobrescu R, Gerber D, Kiss V, Szécsényi-Nagy A, Mende BG, Gallina Z, Somogyi K, Kulcsár G, Gál E, Bendrey R, Allentoft ME, Sirbu G, Dergachev V, Shephard H, Tomadini N, Grouard S, Kasparov A, Basilyan AE, Anisimov MA, Nikolskiy PA, Pavlova EY, Pitulko V, Brem G, Wallner B, Schwall C, Keller M, Kitagawa K, Bessudnov AN, Bessudnov A, Taylor W, Magail J, Gantulga J-O, Bayarsaikhan J, Erdenebaatar D, Tabaldiev K, Mijiddorj E, Boldgiv B, Tsagaan T, Pruvost M, Olsen S, Makarewicz CA, Valenzuela Lamas S, Albizuri Canadell S, Nieto Espinet A, Iborra MP, Lira Garrido J, Rodríguez González E, Celestino S, Olària C, Arsuaga JL, Kotova N, Pryor A, Crabtree P, Zhumatayev R, Toleubaev A, Morgunova NL, Kuznetsova T, Lordkipanize D, Marzullo M, Prato O, Bagnasco Gianni G, Tecchiati U, Clavel B, Lepetz S, Davoudi H, Mashkour M, Berezina NYa, Stockhammer PW, Krause J, Haak W, Morales-Muñiz A, Benecke N, Hofreiter
M, Ludwig A, Graphodatsky AS, Peters J, Kiryushin KYu, Iderkhangai T-O, Bokovenko NA, Vasiliev SK, Seregin NN, Chugunov KV, Plasteeva NA, Baryshnikov GF, Petrova E, Sablin M, Ananyevskaya E, Logvin A, Shevnina I, Logvin V, Kalieva S, Loman V, Kukushkin I, Merz I, Merz V, Sakenov S, Varfolomeyev V, Usmanova E, Zaibert V, Arbuckle B, Belinskiy AB, Kalmykov A, Reinhold S, Hansen S, Yudin AI, Vybornov AA, Epimakhov A, Berezina NS, Roslyakova N, Kosintsev PA, Kuznetsov PF, Anthony D, Kroonen GJ, Kristiansen K, Wincker P, Outram A, Orlando L (2021) The origins and spread of domestic horses from the Western Eurasian steppes. *Nature*, **598**, 634–640. https://doi.org/10.1038/s41586-021-04018-9

Lichtenfels JR, Kharchenko VA, Dvojnos GM (2008) Illustrated identification keys to strongylid parasites (Strongylidae: Nematoda) of horses, zebras and asses (Equidae). *Vet Parasitol*, **156**, 4–161. https://doi.org/10.1016/j.vetpar.2008.04.026

Louro M, Kuzmina TA, Bredtmann CM, Diekmann I, de Carvalho LMM, von Samson-Himmelstjerna G, Krücken J (2021) Genetic variability, cryptic species and phylogenetic relationship of six cyathostomin species based on mitochondrial and nuclear sequences. *Scientific Reports*, **11**, 8245. https://doi.org/10.1038/s41598-021-87500-8

Love MI, Huber W, Anders S (2014) Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. *Genome Biol*, **15**, 550. https://doi.org/10.1186/s13059-014-0550-8

Luu K, Bazin E, Blum MGB (2017) *pcadapt*: an R package to perform genome scans for selection based on principal component analysis. *Molecular Ecology Resources*, **17**, 67–77. https://doi.org/10.1111/1755-0998.12592

Lyons E (2003) Population-S benzimidazole- and tetrahydropyrimidine-resistant small strongyles in a pony herd in Kentucky (1977-1999): effects of anthelmintic treatment on the parasites as determined in critical tests. *Parasitology Research*, **91**, 407–411. https://doi.org/10.1007/s00436-003-0983-6

Lyons ET, Drudge JH, Tolliver SC (1992) Review of prevalence surveys of internal parasites recovered (1951-1990) from horses at necropsy in Kentucky. *Journal of Equine Veterinary Science*, **12**, 9–16. https://doi.org/10.1016/S0737-0806(06)81377-8

Lyons ET, Tolliver SC, Drudge JH (1999) Historical perspective of cyathostomes: prevalence, treatment and control programs. *Vet Parasitol*, **85**, 97–111; discussion 111-2, 215–25. https://doi.org/10.1016/s0304-4017(99)00091-6

Malsa J, Courtot É, Boisseau M, Dumont B, Gombault P, Kuzmina TA, Basiaga M, Lluch J, Annonay G, Dhorne-Pollet S, Mach N, Sutra J-F, Wimel L, Dubois C, Guégnard F, Serreau D, Lespine A, Sallé G, Fleurance G (2022) Effect of sainfoin (*Onobrychis viciifolia*) on cyathostomin eggs excretion, larval development, larval community structure and efficacy of ivermectin treatment in horses. *Parasitology*, **149**, 1439–1449. https://doi.org/10.1017/S0031182022000853

Marini F, Ludt A, Linke J, Strauch K (2021) GeneTonic: an R/Bioconductor package for streamlining the interpretation of RNA-seq data. *BMC Bioinformatics*, **22**, 610. https://doi.org/10.1186/s12859-021-04461-5

Marsh AE, Lakritz J (2023) Reflecting on the past and fast forwarding to present day anthelmintic resistant *Ancylostoma caninum*-A critical issue we neglected to forecast. *International Journal for Parasitology. Drugs and Drug Resistance*, **22**, 36–43. https://doi.org/10.1016/j.ijpddr.2023.04.003

Martin F, Dube F, Karlsson Lindsjö O, Eydal M, Höglund J, Bergström TF, Tydén E (2020) Transcriptional responses in *Parascaris univalens* after in vitro exposure to ivermectin, pyrantel citrate and thiabendazole. *Parasites & Vectors*, **13**, 342. https://doi.org/10.1186/s13071-020-04212-0

Matsubara K, Tarui H, Toriba M, Yamada K, Nishida-Umehara C, Agata K, Matsuda Y (2006) Evidence for different origin of sex chromosomes in snakes, birds, and mammals and step-wise differentiation of snake sex chromosomes. *Proceedings of the National Academy of Sciences*, **103**, 18190–18195. https://doi.org/10.1073/pnas.0605274103
McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, Garimella K, Altshuler D, Gabriel S, Daly M, DePristo MA (2010) The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. *Genome Res*, **20**, 1297–303.
https://doi.org/10.1101/gr.107524.110

Merchant M, Mata CP, Liu Y, Zhai H, Protasio AV, Modis Y (2022) A bioactive phlebovirus-like envelope protein in a hookworm endogenous virus. *Science Advances*, **8**, eabj6894.
https://doi.org/10.1126/sciadv.abj6894

Mistry J, Chuguransky S, Williams L, Qureshi M, Salazar GA, Sonnhammer ELL, Tosatto SCE, Paladin L, Raj S, Richardson LJ, Finn RD, Bateman A (2021) Pfam: The protein families database in 2021. *Nucleic Acids Research*, **49**, D412–D419. https://doi.org/10.1093/nar/gkaa913

Morgulis A, Gertz EM, Schäffer AA, Agarwala R (2006) A fast and symmetric DUST implementation to mask low-complexity DNA sequences. *Journal of Computational Biology: A Journal of Computational Molecular Cell Biology*, **13**, 1028–1040. https://doi.org/10.1089/cmb.2006.13.1028

Moser W, Schindler C, Keiser J (2017) Efficacy of recommended drugs against soil transmitted helminths: systematic review and network meta-analysis. *BMJ (Clinical research ed.)*, **358**, j4307.
https://doi.org/10.1136/bmj.j4307

Neveu C, Charvet CL, Fauvin A, Cortet J, Beech RN, Cabaret J (2010) Genetic diversity of levamisole receptor subunits in parasitic nematode species and abbreviated transcripts associated with resistance. *Pharmacogenet Genomics*, **20**, 414–25.
https://doi.org/10.1097/FPC.0b013e328338ac8c

Nielsen MK (2022) Anthelmintic resistance in equine nematodes: Current status and emerging trends. *International Journal for Parasitology. Drugs and Drug Resistance*, **20**, 76–88.
https://doi.org/10.1016/j.ijpddr.2022.10.005

Nielsen MK, Steuer AE, Anderson HP, Gavriliuc S, Carpenter AB, Redman EM, Gilleard JS, Reinemeyer CR, Poissant J (2022) Shortened egg reappearance periods of equine cyathostomins following ivermectin or moxidectin treatment: morphological and molecular investigation of efficacy and species composition. *International Journal for Parasitology*, **52**, 787–798.
https://doi.org/10.1016/j.ijpara.2022.09.003

Ogbourne CP (1972) Observations on the free-living stages of strongylid nematodes of the horse. *Parasitology*, **64**, 461–477. https://doi.org/10.1017/S0031182000045534

Patro R, Duggal G, Love MI, Irizarry RA, Kingsford C (2017) Salmon provides fast and bias-aware quantification of transcript expression. *Nature Methods*, **14**, 417–419.
https://doi.org/10.1038/nmeth.4197

Patterson N, Moorjani P, Luo Y, Mallick S, Rohland N, Zhan Y, Genschoreck T, Webster T, Reich D (2012) Ancient admixture in human history. *Genetics*, **192**, 1065–93.
https://doi.org/10.1534/genetics.112.145037

Paysan-Lafosse T, Blum M, Chuguransky S, Grego T, Pinto BL, Salazar GA, Bileschi ML, Bork P, Bridge A, Colwell L, Gough J, Haft DH, Letunić I, Marchler-Bauer A, Mi H, Natale DA, Orengo CA, Pandurangan AP, Rivoire C, Sigrist CJA, Sillitoe I, Thanki N, Thomas PD, Tosatto SCE, Wu CH, Bateman A (2023) InterPro in 2022. *Nucleic Acids Research*, **51**, D418–D427.
https://doi.org/10.1093/nar/gkac993

Poissant J, Gavriliuc S, Bellaw J, Redman EM, Avramenko RW, Robinson D, Workentine ML, Shury TK, Jenkins EJ, McLoughlin PD, Nielsen MK, Gilleard JS (2021) A repeatable and quantitative DNA metabarcoding assay to characterize mixed strongyle infections in horses. *International Journal for Parasitology*, **51**, 183–192. https://doi.org/10.1016/j.ijpara.2020.09.003

Pollet N (2025) Genomic and transcriptomic insights into the genetic basis of anthelmintic resistance in a cyathostomin parasitic nematode. *Peer Community in Genomics*, 100290.
https://doi.org/10.24072/pci.genomics.100290
Raineri E, Ferretti L, Esteve-Codina A, Nevado B, Heath S, Pérez-Enciso M (2012) SNP calling by sequencing pooled samples. *BMC Bioinformatics*, **13**, 239. https://doi.org/10.1186/1471-2105-13-239

Robertson AP, Bjørn HE, Martin RJ (2000) Pyrantel resistance alters nematode nicotinic acetylcholine receptor single-channel properties. *European Journal of Pharmacology*, **394**, 1-8.
https://doi.org/10.1016/s0014-2999(00)00135-7

Robertson SJ, Pennington AJ, Evans AM, Martin RJ (1994) The action of pyrantel as an agonist and an open channel blocker at acetylcholine receptors in isolated Ascaris suum muscle vesicles. *Eur J Pharmacol*, **271**, 273-82. https://doi.org/10.1016/0014-2999(94)90784-6

Sallé G, Cortet J, Bois I, Dubes C, Guyot-Sionest Q, Larrieu C, Landrin V, Majorel G, Wittreck S, Woringer E, Courouce A, Guillot J, Jacquiet P, Guegnard F, Blanchard A, Leblond A (2017) Risk factor analysis of equine strongyle resistance to anthelmintics. *Int J Parasitol Drugs Drug Resist*, **7**, 407-415. https://doi.org/10.1016/j.ijpddr.2017.10.007

Sallé G, Doyle SR, Cortet J, Cabaret J, Berriman M, Holroyd N, Cotton JA (2019) The global diversity of *Haemonchus contortus* is shaped by human intervention and climate. *Nature Communications*, **10**, 4811. https://doi.org/10.1038/s41467-019-12695-4

Sallé G, Guillot J, Tapprest J, Foucher N, Sevin C, Laugier C (2020) Compilation of 29 years of postmortem examinations identifies major shifts in equine parasite prevalence from 2000 onwards. *International Journal for Parasitology*, **50**, 125-132. https://doi.org/10.1016/j.ijpara.2019.11.004

Sallé G, Kornaś S, Basiaga M (2018) Equine strongyle communities are constrained by horse sex and species dipersal-fecundity trade-off. *Parasites & Vectors*, **11**, 279. https://doi.org/10.1186/s13071-018-2858-9

von Samson-Himmelstjerna G, Thompson RA, Krücken J, Grant W, Bowman DD, Schnyder M, Deplazes P (2021) Spread of anthelmintic resistance in intestinal helminths of dogs and cats is currently less pronounced than in ruminants and horses – Yet it is of major concern. *International Journal for Parasitology: Drugs and Drug Resistance*, **17**, 36-45.
https://doi.org/10.1016/j.ijpddr.2021.07.003

Sargison N, Chambers A, Chaudhry U, Costa Júnior L, Doyle SR, Ehimiyein A, Evans M, Jennings A, Kelly R, Sargison F, Sinclair M, Zahid O (2022) Faecal egg counts and nemabiome metabarcoding highlight the genomic complexity of equine cyathostomin communities and provide insight into their dynamics in a Scottish native pony herd. *International Journal for Parasitology*, **52**, 763-774.
https://doi.org/10.1016/j.ijpara.2022.08.002

Scare JA, Lyons ET, Wielgus KM, Nielsen MK (2018) Combination deworming for the control of double-resistant cyathostomin parasites - short and long term consequences. *Veterinary Parasitology*, **251**, 112-118. https://doi.org/10.1016/j.vetpar.2018.01.010

Schwarz EM, Hu Y, Antoshechkin I, Miller MM, Sternberg PW, Aroian RV (2015) The genome and transcriptome of the zoonotic hookworm *Ancylostoma ceylanicum* identify infection-specific gene families. *Nature Genetics*, **47**, 416-422. https://doi.org/10.1038/ng.3237

Seppey M, Manni M, Zdobnov EM (2019) BUSCO: Assessing Genome Assembly and Annotation Completeness. *Methods in Molecular Biology* (Clifton, N.J.), **1962**, 227-245.
https://doi.org/10.1007/978-1-4939-9173-0_14

Sieburth D, Ch'ng Q, Dybbs M, Tavazoie M, Kennedy S, Wang D, Dupuy D, Rual J-F, Hill DE, Vidal M, Ruvkun G, Kaplan JM (2005) Systematic analysis of genes required for synapse structure and function. *Nature*, **436**, 510-517. https://doi.org/10.1038/nature03809

da Silva Anjos DH, de Lurdes A Rodrigues M (2006) Diversity of the infracommunities of strongylid nematodes in the ventral colon of *Equus caballus* from Rio de Janeiro state, Brazil. *Vet Parasitol*, **136**, 251-257. https://doi.org/10.1016/j.vetpar.2005.05.070
Silva AV, Costa HM, Santos HA, Carvalho RO (1999) Cyathostominae (Nematoda) parasites of Equus caballus in some Brazilian states. *Vet Parasitol*, **86**, 15-21. https://doi.org/10.1016/s0304-4017(99)00078-3

Sleigh JN (2010) Functional analysis of nematode nicotinic receptors. *Bioscience Horizons*, **3**, 29-39.
https://doi.org/10.1093/biohorizons/hzq005

Tarailo-Graovac M, Chen N (2009) Using RepeatMasker to identify repetitive elements in genomic sequences. *Current Protocols in Bioinformatics*, **Chapter 4**, 4.10.1-4.10.14.
https://doi.org/10.1002/0471250953.bi0410s25

The cyathostomin genomics consortium (2023) Cyathostomin genomics. https://doi.org/10.57745/1FB3SD

Tombak KJ, Hansen CB, Kinsella JM, Pansu J, Pringle RM, Rubenstein DI (2021) The gastrointestinal nematodes of plains and Grevy's zebras: Phylogenetic relationships and host specificity. *International Journal for Parasitology: Parasites and Wildlife*, **16**, 228-235.
https://doi.org/10.1016/j.ijppaw.2021.10.007

Tydén E, Enemark HL, Franko MA, Höglund J, Osterman-Lind E (2019) Prevalence of *Strongylus vulgaris* in horses after ten years of prescription usage of anthelmintics in Sweden. *Veterinary Parasitology*: *X*, **2**, 100013. https://doi.org/10.1016/j.vpoa.2019.100013

Uliano-Silva M, Ferreira JGRN, Krasheninnikova K, Darwin Tree of Life Consortium, Formenti G, Abueg L, Torrance J, Myers EW, Durbin R, Blaxter M, McCarthy SA (2023) MitoHiFi: a python pipeline for mitochondrial genome assembly from PacBio high fidelity reads. *BMC bioinformatics*, **24**, 288. https://doi.org/10.1186/s12859-023-05385-y

Wang C, Torgerson PR, Kaplan RM, George MM, Furrer R (2018) Modelling anthelmintic resistance by extending eggCounts package to allow individual efficacy. *International Journal for Parasitology. Drugs and Drug Resistance*, **8**, 386-393. https://doi.org/10.1016/j.ijpddr.2018.07.003

Welz C, Kruger N, Schniederjans M, Miltsch SM, Krucken J, Guest M, Holden-Dye L, Harder A, von Samson-Himmelstjerna G (2011) SLO-1-channels of parasitic nematodes reconstitute locomotor behaviour and emodepside sensitivity in *Caenorhabditis elegans* slo-1 loss of function mutants. *PLoS Pathog*, **7**, e1001330. https://doi.org/10.1371/journal.ppat.1001330

Wingett SW, Andrews S (2018) FastQ Screen: A tool for multi-genome mapping and quality control. *F1000Research*, **7**, 1338. https://doi.org/10.12688/f1000research.15931.2

Zhang Y, Park C, Bennett C, Thornton M, Kim D (2021) Rapid and accurate alignment of nucleotide conversion sequencing reads with HISAT-3N. *Genome Research*, **31**, 1290-1295.
https://doi.org/10.1101/gr.275193.120