The identification of small molecule inhibitors with anthelmintic activities that target conserved proteins among ruminant gastrointestinal nematodes

ABSTRACT Gastrointestinal nematode (GIN) infections are a major concern for the ruminant industry worldwide and result in significant production losses. Naturally occurring polyparasitism and increasing drug resistance that potentiate disease outcomes are observed among the most prevalent GINs of veterinary importance. Within the five major taxonomic clades, clade Va represents a group of GINs that predominantly affect the abomasum or small intestine of ruminants. However, the development of effective broad-spectrum anthelmintics against ruminant clade Va GINs has been challenged by a lack of comprehensive druggable genome resources. Here, we first assembled draft genomes for three clade Va species (Cooperia oncophora, Trichostrongylus colubriformis, and Ostertagia ostertagi) and compared them with closely related ruminant GINs. Genome-wide phylogenetic reconstruction showed a relationship among ruminant GINs structured by taxonomic classification. Orthogroup (OG) inference and functional enrichment analyses identified 220 clade Va-specific and Va-conserved OGs, enriched for functions related to cell cycle and cellular senescence. Further transcriptomic analysis identified 61 taxonomically and functionally conserved clade Va OGs that may function as drug targets for new broad-spectrum anthelmintics. Chemogenomic screening identified 11 compounds targeting homologs of these OGs, thus having potential anthelmintic activity. In in vitro phenotypic assays, three kinase inhibitors (digitoxigenin, K-252a, and staurosporine) exhibited broad-spectrum anthelmintic activities against clade Va GINs by obstructing the motility of exsheathed L3 (xL3) or molting of xL3 to L4. These results demonstrate valuable applications of the new ruminant GIN genomes in gaining better insights into their life cycles and offer a contemporary approach to discovering the next generation of anthelmintics. IMPORTANCE Gastrointestinal nematode (GIN) infections in ruminants are caused by parasites that inhibit normal function in the digestive tract of cattle, sheep, and goats, thereby causing morbidity and mortality. Coinfection and increasing drug resistance to current therapeutic agents will continue to worsen disease outcomes and impose significant production losses on domestic livestock producers worldwide. In combination with ongoing therapeutic efforts, advancing the discovery of new drugs with novel modes of action is critical for better controlling GIN infections. The significance of this study is in assembling and characterizing new GIN genomes of Cooperia oncophora, Ostertagia ostertagi, and Trichostrongylus colubriformis for facilitating a multi-omics approach to identify novel, biologically conserved drug targets for five major GINs of veterinary importance. With this information, we were then able to demonstrate the potential of commercially available compounds as new anthelmintics.

G astrointestinal nematode (GIN) infections are one of the major problems facing livestock producers worldwide, and there is growing concern about the increase in drug resistance to all major classes of anthelmintics (1,2).Once controlled predomi nantly by anthelmintics, overuse and abuse have led to resistance against most major classes of drugs (3)(4)(5)(6)(7).Meat and milk producers have resorted to less than adequate management-based control strategies, especially for those parasites that are the most pathogenic, i.e., Haemonchus spp., Ostertagia ostertagi, and Teladorsagia circumcincta.Although drug resistance in Ostertagia is in its early stages, biological similarities to Teladorsagia suggest that greater dissemination of Ostertagia-resistant populations is not far off.The discovery and development of new effective drugs have been slow; consequently, some regions of the world have adopted techniques involving refugia along with more disciplined drug intervention to control diseases caused by the most prolific and pathogenic nematodes (8)(9)(10)(11)(12)(13)(14).In other locations, producers and manufactur ers have been working to extend the therapeutic lifespan of the current drugs by using combination or rotational therapies, though in the US, combination therapy has been hindered by government policy and regulation (15)(16)(17)(18).
In trichostrongylid nematodes, drug resistance is likely derived not from the introduction of new genetic mutations but instead by selecting pre-existing mutations in genetically diverse worm populations (19).The rise in resistant worms in drug-treated hosts depends upon the frequency and level of drug intervention, combined with opportunity, conditions, the allele frequency of resistant SNPs, and "sloppy fitness" [the evolution of phenotypic plasticity and flexibility to rapidly adapt to novel condi tions and/or hosts (20)].Over time, the selection of a pre-existing subpopulation of drug-resistant worms leads to a shift in the overall population in favor of the resistant genotype (21).
In general, each class of anthelmintics has enjoyed a life span of 15-20 years before resistance has surfaced (22).This is related to the degree of complexity in the resistance pathway and the large effective worm population size that leads to high intra-species and intra-population genetic variability.In addition, the problem is exacerbated by extensive gene flow among populations and/or the plasticity of the target gene (23).For example, β-tubulin (the target for benzimidazoles) is highly conserved among organisms, whereas the high variability in nicotinic acetylcholine receptors (e.g., acr-23), the binding moiety for amino-acetonitrile derivatives (Monepantel), likely hastens the development of resistance to these drugs (24).Thus, if drug treatment is to continue to be the mode of GIN control going forward, chemicals with novel modes of action are needed in combination with integrated pasture management.
Genome and transcriptome sequencing have demonstrably increased in recent years, including those of larger genomes, thanks to state-of-the-art, next-generation sequenc ing techniques.Searching for genetic commonalities among parasitic helminths has, to date, adopted the paradigm "the more the merrier" by comparing as many of the available databases to identify orthologous or pan-nematode target sequences (25).The International Helminth Genome Consortium (25) evaluated draft genomes from 81 parasitic and non-parasitic worms and identified links to parasitism and new drug targets among highly diverse organisms.Their comparisons characterized far-reaching, lineage-specific differences in core metabolism-related genes and protein families, some of which had previously been targeted for drug development.Using an in silico screen, the authors also identified new promising drug targets from parasitic nematodes and platyhelminths and exploited a drug repurposing approach to prioritize potential anthelmintics.Recently, 409 of 817 proposed drugs were experimentally screened against the parasitic nematode Trichuris muris, and 50 of 409 drugs (~12%) showed complete or partial motility inhibition against adult worms (26).However, these efforts have not yielded candidates for broad-spectrum anthelmintics; sequence homology and orthology among disparate organisms can be complicated by exaptation and co-option (the repurposing of genes/proteins with similar sequences to serve different functions among dissimilar pathogens).Thus, the activities of genes among evolutionarily distinct lineages may not be congruent, which in turn poses challenges for identifying pan-nem atode targets for treatment.
To circumvent potential issues imposed by evolution and genetic diversity and identify drug or vaccine targets shared by ruminant GIN, we chose to cull our data base to the genomes of a select group of Trichostrongyloidea that more closely share an evolutionary history, are commonly found among livestock hosts, and pose an economic challenge to producers worldwide.To facilitate this approach, we first sequenced, assembled, and annotated three Trichostrongyloidea genomes of major veterinary importance (two small intestinal parasites: Cooperia oncophora, Trichostrongy lus colubriformis, and one abomasal parasite: O. ostertagi) and compared them with two existing clade Va ruminant abomasal parasites [Haemonchus contortus (27) and Teladorsagia circumcincta (28)].Other clade V members belonging to the Chabertii dae [the large intestinal parasite: Oesophagostomum dentatum in clade Vc (29)] and Dictyocaulidae [the lung parasite: Dictyocaulus viviparus in clade Vb (30)] were selec ted as outgroups.Our primary goal was identifying biologically conserved targets in the selective group of small intestinal/abomasal parasites in ruminants (clade Va) for therapeutic intervention.We hypothesized that conservation in the evolution of parasitism is constrained by stratified ecological niches within the hosts, with the small intestine/abomasum serving for digestion and nutrient absorption, as opposed to other compartments like the large intestine and lung.We performed comparative genomics and transcriptomic analyses of our interest group with the outgroup.This enabled the identification of gene families specific to and/or shared among ruminant GINs (small intestine and abomasum) in clade Va and provided a better understanding of transcriptional conservation within the selected gene families during their complex life cycles.These findings were extended beyond identifying conserved and functional drug target candidates to performing chemogenomic screening using the ChEMBL database to prioritize commercially available compounds.We then evaluated existing drugs as potential anthelmintics with broad-spectrum activity against ruminant GINs in clade Va and identified three small molecular inhibitors exhibiting promising drug activities.This approach shows promise to leverage comparative genomics for drug discovery and translational medicine.

The genomes of C. oncophora, T. colubriformis, and O. ostertagi reveal genomic features and infer phylogenetic relationship driven by taxonomic classification
To facilitate omics-driven drug discovery approaches on the major GIN species of importance to cattle, sheep, and goats, we undertook a stepwise approach (Fig. 1).First, we sequenced and generated draft genomes of three Trichostrongyloidea species (C.oncophora, T. colubriformis, and O. ostertagi) with no existing genomic resources.This greatly expanded the genetic diversity of the clade Va species for comparative genomic analyses.The draft genomes were then annotated per gene using KEGG, Gene Ontology (GO), and InterPro databases and compared with previously generated genomes of other ruminant parasites in clade V (Table 1).Following removal of redundant contigs, the sizes of the three newly assembled genomes were estimated to be 476.The assembly metrics including the number of contigs, L50, and N50 statistics indicated that T. colubriformis was the most contiguous assembly among the three newly sequenced species and comparable to T. circumcincta and O. dentatum.The genome of C. oncophora had by far the highest number of protein-coding genes (n = 22,571) among the seven clade V species, but its gene density (47.4 genes/Mb) was lower than the chromosome-level genome of H. contortus (68.7 genes/Mb) and the smallest but most compact genome of D. viviparus (81.9 genes/Mb).CEGMA-based completeness estimated by 248 highly conserved core eukaryote gene sets, ranged between 86.7% and 91.5%.The BUSCO completeness estimated from 3,131 nematode markers was between 76.1% and 83.2% when partial matches were included.Consistent with the assembly metrics, T. colubriformis was the next most complete genome assembly following H. contortus and D. viviparus based on CEGMA score (91.5%).The transcriptome mapping rates, including both unique and multiple mapped RNA sequencing (RNA-seq) reads, were as follows: C. oncophora (84.3%), T. colubriformis (82.7%),O. ostertagi (70.4%), H. contortus (93.4%), and T. circumcincta (84.9%) (Table S3).These observations validated that the new draft genomes are of similar quality to previously published genomes of other clade V species and are therefore suitable for comparative genomic analyses (Fig. S1).All genome assemblies and genomic and transcriptomic sequencing data sets generated and/or used in this study are summarized in Tables S1 to S3.
Next, we performed whole-genome-based clustering that specifically inferred phylogenetic relationships among the major ruminant parasites in clade V (Fig. 2).A species tree was generated by the STAG algorithm and rooted by the STRIDE algorithm used in OrthoFinder, a tool for phylogenetic orthology inference (31).The genome-scale phylogenetic tree, based on 5,181 orthologous groups comprising genes from our subset of clade V nematodes, revealed three distinct clades, clade Va, Vb, and Vc as previously observed, placing O. dentatum, the large intestine parasite (clade Vc) in the earliest branching clade (25).The group of five clade Va species residing in the abomasum or small intestine was divided into two subgroups.These relationships corroborate previous phylogenetic inferences derived from studies comparing mito chondrial genomes, a subset of protein-coding genes, and morphological characteristics of Trichostrongyloidea; O. dentatum (nodular worm in large intestine in Strongyloidea superfamily) and D. viviparus (lungworm in Trichostrongyloidea superfamily) comprised a distinct clade separated from other GINs in trichostrongyloid nematodes (32)(33)(34).

Orthogroups in clade Va species show taxonomic and functional conserva tion and diversification
The taxonomic distribution of genes across ruminant GINs was investigated to better understand the level of conservation and diversification of clade Va orthologous groups.We used this analysis to prioritize clade Va conserved genes (proteins) necessary for parasite survival, thereby broadly vulnerable to drugs against ruminant GIN infections.This is extremely important considering the prevalence and polyparasitism of these trichostrongylids in ruminant hosts (35)(36)(37)(38) and the increasing propensity of resistance to currently available anthelmintics.Homology relationships between protein-coding genes spanning trichostrongylid nematodes were inferred using an OG inference algorithm (31).An OG refers to a set of genes descended from a single ancestral gene of a given set of species, which includes both orthologs and paralogs.The OGs in trichostrongylids were defined using the five ruminant GINs (clade Va: C. oncophora, T. colubriformis, O. ostertagi, H. contortus, and T. circumcincta) residing in the abomasum or small intestine as taxonomic groups of interest, with D. viviparus in the lung (clade Vb) and O. dentatum in the large intestine (clade Vc) serving as outgroup species (Fig. 2).About 90% of genes from each species shared homology among the clade V species, comprising a total of 17,689 OGs (Fig. 3A).The remaining 10% were species-specific "singletons" unassigned to any OG.After predicting functional annotations of genes in each species, putative functions of OGs were determined if 50% or more gene members of that OG agreed with their func tional annotation(s) (Table S4).Functional enrichment analysis was performed with the consensus annotations of OGs.
A total of 6,486 OGs were conserved in all examined clade Va species, of which 79.9% (5,181 OGs) were also shared with the outgroup species, suggesting fundamental roles for these genes in the ruminant parasites in clade V (Fig. 3A; Fig. S2).Functional enrichment analysis among the clade Va-conserved (VaC) OGs (relative to all OGs; Fig. 3B) identified enrichment of transcription factors, mRNA/ribosome/mitochondria biogenesis, kinases and signaling pathways involved in cell growth, differentiation, apoptosis, and homeostasis.These functions mediate various cellular processes and are often found in not only parasitic nematodes but also all kingdoms of living organisms.
A total of 4,706 clade Va-specific (VaS) OGs were identified to have OG membership supported by one or more clade Va species but not by the outgroup species.These OGs represented an average of 15.3%-18% of genes from each species, except for H. contortus (with 23.6%), an observation suggesting missing genes in the non-Haemonchus genomes due either to technical limitations with the assembly/annotations or to true biological loss of the genes (Fig. 3A).The most significantly enriched term among VaS OGs was a pathway related to proteases and protease inhibitors (Fig. 3C).This may link their roles to unique parasitic life styles involving immunomodulation, host tissue penetration, and nutrient digestion, making them favorable drug target candidates (39,40).Proteases were one of the protein families that displayed frequent gene gains and losses among the GINs in clade Va, where metallopeptidases appeared to be a rapidly evolving protease family not only in the clade Va species but also in the outgroup species (Fig. S3).The expansion of metallopeptidases was previously reported in clade IVa parasites whose invasions necessitate skin penetration and/or degradation of tissues while migrating through the digestive tracts and other organs (25,41).On the other hand, differential expansion of aspartic proteases was quite exclusive to the clade Va species, suggesting their specialized roles in GINs in ruminant hosts; however, the draft nature of the genomes cannot be ruled out as a contributing factor (Fig. S3).One of the interesting InterPro protein domains enriched in the VaS OGs was "secreted nematode clade V proteins" (IPR035126; FDR ≤ 1.04 × 10 −5 ), which was greatly expanded in GIN species (clades Va and Vc) but not in the lung parasite, D. viviparus (clade Vb).In our data set, 19 OGs corresponding to 206 orthologous genes were annotated with "secreted clade V protein" by InterPro and most of them were predicted to have signal peptides by either SignalP (42) or SecretomeP (43) (Table S4).The secreted clade V domain was previously shown to be abundant in Strongylida parasites, and their transcriptional expressions were significantly upregulated after infection (44).Thus, this protein family may have evolved to aid parasitism in the hosts during the infective stages and therefore represents a new class of parasite-specific targets for drug development.
A total of 220 clade Va-specific and Va-conserved (VaSC) OGs (1.24% of total OGs) were identified based on gene members being present in all five clade Va species but absent in the two outgroup species (Fig. 3A).Functional enrichment analysis of the 220 VaSC OGs identified pathways related to cell cycle, apoptosis, the p53 signaling pathway, and protein glycosylation, which partially overlapped with the enriched pathways in either the conserved or unique OGs of clade Va species (Fig. 3D).The most significantly enriched pathway was "cellular senescence" (map04218; P ≤ 0.0013), related to arrest of cellular development and proliferation, which is likely attributed to shared evolutionary constraints such as dynamic and hostile environments that all clade Va species encoun ter during two separate phases of their life: free-living and parasitic stages.Of note, approximately 18% of the 220 VaSC OGs was single-copy orthologs and 60% of them comprised a set of just one or two copies of genes from each species (Table S4).This suggests that most of the clade Va exclusively conserved OGs remained as one or two copies, with minimal duplication or loss events since the last common ancestor.The hypothesis that functional conservation among OGs can be used to identify promising drug targets was further supported by results from homology inference with three ruminant host species (cattle, sheep, and goat), showing that 182 of the total 220 VaSC OGs (82.7%) were missing from all of these hosts (Table S5).Taken together, 220 VaSC OGs seem to be ideal candidates as anthelmintic drug targets and therefore were selected for further prioritization.

Stage-specific comparative transcriptomic analysis determined interspecies conservation and essential functionality of 220 VaSC OGs in clade Va
In addition to the taxonomic-and sequence homology-based approach, we analyzed the stage-specific gene expression profiles of the 220 VaSC OGs (2,221 genes), expanding our understanding of the transcriptional regulation of genes that have likely evolved in the trichostrongylid species.For comparative transcriptomic analysis of the 220 VaSC OGs, RNA sequencing data sets were collected from publicly or internally produced sources, considering expression profiles from both free-living and parasitic stages (Table S3).To ascertain expression within a species, each gene was evaluated for constitutive (CEG) or differential (DEG) expression (Fig. 4).If a gene had fragments per kilobase per million reads (FPKM) ≥ 1 in (i) all available, (ii) only free-living, or (iii) only parasitic stages in a corresponding species, it was annotated as a "CEG in all, " a "CEG in free-living, " or a "CEG in parasitic, " respectively (Fig. 5A).If the differential expression of a gene was statistically significant between free-living and parasitic stages, the gene was designated as either a "higher DEG in free-living" or a "higher DEG in parasitic" (Fig. 5B).These data were used to target parasite genes that are potentially essential in parasitic stages, so DEGs with significant overexpression in parasitic stages compared to free-living stages ("higher DEGs in parasitic stages") meet the criteria for consideration as drug target candidates.In addition to DEGs, CEGs, whether or not they overlap with DEGs, could also be regarded as potential drug candidates due to their predicted importance in maintaining essential cellular functions in parasitic stages.For instance, genes classified as "CEGs in parasites" are genes that exhibit continuous expression throughout the parasitic stages, suggesting critical roles for survival within the host.By categorizing these genes according to our defined criteria, we aim to evaluate a comprehensive set of genes, including both CEGs and DEGs, which could be explored as potential druggable targets.Gene expression within a species was summarized for each OG (Fig. 4).Among the 220 VaSC OGs, each consisted of one or more gene members that originated from each of the five clade Va species.Conserved expression at the OG level was determined based on the presence of at least one or more representative genes being expressed in each of the clade Va species.If the representative genes from each of the five clade Va species agreed on one of the five expression categories described above, the OG was defined as having conserved cross-species gene expression and was assigned to the corresponding category: constitutively expressed OG in all (CEO in all), free-living (CEO in free-living), or parasitic (CEO in parasitic) stages, or differentially expressed OG, higher in free-living (higher DEO in free-living), or parasitic (higher DEO in parasitic) stages (Fig. 4).
The highest proportion of genes and OGs were classified as CEGs and CEOs in all stages, respectively, followed by CEGs/CEOs in free-living, CEGs/CEOs in parasitic, higher DEGs/DEOs in parasitic, and higher DEGs/DEOs in free-living stages (Fig. 5).Among the 220 VaSC OGs, 206 OGs (93.6%) presented evidence of transcription at either free-living or parasitic stages from at least one of the five species, and 61 OGs (27.7%) showed conserved expression patterns such as CEOs and DEOs across all clade Va species (Fig. 5C; Table S6).Of the 58 CEOs, 39 showed CEOs in all, 13 were CEOs in free-living, and 6 were CEOs in parasitic stages.One of the CEOs in all, OG0000036, had the second-largest gene membership among the 220 VaSC OGs; 45 of the 60 genes were expressed in both free-living and parasitic stages (Table S6).However, its putative function could not be determined due to a lack of annotated genes in the group; functions of the genes even in best annotated species, H. contortus, were not able to be inferred from available protein databases.Given that the genes comprising the 220 VaSC OGs are conserved among clade Va only, and without homology with the ruminant hosts, we putatively defined them as parasite-specific genes whose suitability as drug targets deserves evaluation.Intriguingly, there were no OGs with significantly higher expression in the free-living stage (higher DEO in free-living), but three OGs were more highly expressed in the parasitic stage (higher DEO in parasitic) (Fig. 5C).Of the three DEOs higher in the parasitic stage, two had consensus functional anno tations substantiated by more than 50% of gene members in each OG: OG0002618 as "globin-like superfamily" and OG0009146 as "alpha crystallin/heat-shock protein 20 (hsp20)." HSP20 is a heat shock protein family of molecular chaperones playing various functional roles including response to sudden changes in environment.The α-crystallin domain is conserved among most organisms, but the overall sequence of hsp20 genes remains highly divergent (45,46).A dendrogram of all hsp20 genes identified in the seven clade V species and non-parasitic nematode, C. elegans, resulted in two distinct groups: one specific to the parasites and a second that included C. elegans (Fig. S4).Overall, the transcriptional analysis allowed us to estimate the degree of expression of the trichostrongylid candidate genes that can be utilized as surrogates for their protein expression.In addition, these data revealed that more than one-quarter of the 220 VaSC OGs exhibited interspecies conservation of gene expression in free-living and/or parasitic stages, suggesting functional dependency.This helped further narrow the drug target candidates to those having broad-spectrum anthelmintic potential for GIN infections.

Chemogenomic, data-driven prioritization, identification, and repurposing of drugs with anthelmintic potential
One approach to identifying and prioritizing compounds with anthelmintic potential is to repurpose drugs that are commercially and readily available and previously tested in other species.This approach provides considerable information for predicting efficacy and safety against new parasite targets.With a focus on the 220 VaSC OGs, we performed chemogenomic screening to associate the trichostrongylid genes (proteins) to drug targets in the ChEMBL database predicated on sequence similarity and to prioritize drug-like compounds having potential broad-spectrum activity against the GINs (47) (Fig. 6).First, a BLAST search was performed using the 2,221 protein-coding genes against 9,211 ChEMBL targets with an e-value cutoff ≤ 10 −4 ; within each OG, an agreement of target hits across the five clade Va species was evaluated.This process yielded 882 consensus ChEMBL targets including transcription factors, chaperones, glycosyltransfera ses, dehydrogenases, kinases, transporters, nematode cuticle collagens, and proteases with the potential for broad control or specific application as anthelmintics.
Following the selection of homologous targets in the ChEMBL database, target-com pound pair information and affinity scores (pChEMBL scores ≥ 5) were taken into consideration (Fig. 6).The stepwise prioritization identified 82,812 target-compound pairs, presumably related to the conserved parasite target genes spanning 26 OGs.Combining the expression patterns of 220 VaSC OGs (CEOs or DEOs) obtained from the previous step led to the selection of 10 OGs and associated 57,460 candidates as repositioned anthelmintics.Furthermore, the compounds were ranked based on weighted quantitative estimates of drug likeness (QED) and pChEMBL scores; up to 20 potential drug-like compounds were further culled within each OG resulting in 143 candidates.By determining which of these compounds were commercially available guided by the ZINC database (48), we were able to finalize 11 hit compounds spanning nine OGs with high potential activity against parasite gene products (Table 2).It is worth noting that the ChEMBL target proteins of the final candidate compounds correlated well with the predicted functional roles of the nine OGs, which include kinases, proteases, a heat shock protein, deubiquitinase, dehydrogenase, and a nucleotide-binding protein (Table S7).

Experimental screening of the prioritized compounds against exsheathed L3 of clade Va species
The 11 candidate compounds were screened in vitro against C. oncophora, T. colubrifor mis, and O. ostertagi exsheathed L3 (xL3) mimicking the first parasitic stage following ingestion of the free-living sheathed L3 (L3) by their hosts.First, a single-dose treatment at 100 µM was performed to estimate the efficacy of the 11 compounds based on relative motility inhibition of drug-treated worms; 1% solvent (either DMSO or methanol)-treated worms were used as controls (Table S8).Phenotypic screening showed that 3 of the 11 prioritized compounds induced discernible decreases in movements of xL3 of all three species by day 3 post-treatment, consistent with the predicted broad-spectrum activity on GINs (Fig. 7A through C; Fig. S5).Overall, the intestinal parasites, C. oncophora and T. colubriformis, were more susceptible to the drug treatments than the abomasum parasite, O. ostertagi.The three active compounds digitoxigenin, K-252a, and staurospor ine caused significant decreases (≥75% inhibition) in the motility of C. oncophora and T. colubriformis xL3, and moderate reduction (≥50%) in the motility against O. ostertagi xL3 at a concentration of 100 µM.
Next, IC 50 values of the three active compounds were investigated to determine the potency of the drugs.K-252a, a serine/threonine kinase inhibitor, was the most potent; IC 50 values were 4.4 µM in both C. oncophora and T. colubriformis and 20.2 µM in O. ostertagi (Fig. 7D).Consistent with the previous observation, estimated IC 50 values of the three compounds were the lowest in T. colubriformis, which was the most vulnerable to the drug treatments, followed by C. oncophora and O. ostertagi.We also evaluated the relative motility inhibition of H. contortus xL3 at 100 µM as a readout to estimate the drug efficacy.Although the motility of the drug-treated worms was not significantly different from the motility of the control (data not shown), the phenotype and developmental progress of the H. contortus worms were noticeably affected (Fig. S6).After 24 hours of exposure to staurosporine, the xL3 exhibited an "outstretched" phenotype, and a significant proportion of these larvae had this distinct phenotype on day 3 (Fig. S6J) compared to controls, which were populated with S-or curly shaped larvae (Fig. S6N).In addition, staurosporine-and K-252a-treated H. contortus xL3 completely failed to molt to L4 by day 7.In the digitoxigenin-treated worms and the DMSO control, we observed empty cuticles presumably shed from xL3, suggesting that molting to L4 had occurred by day 7 (indicated by red arrows in Fig. S6C and O); empty cuticles were not observed in K-252a-and staurosporine-treated worms on day 7 (Fig. S6G and S5K).The outstretched phenotype and the molting defect in H. contortus after the staurosporine treatment have not previously been described, even though the motility inhibitory effects of staurosporine and K-252a were previously reported in L1, xL3, and/or L4 H. contortus (49)(50)(51).Of note, this could be attributed to the differences in experimental design (the number of worms per replicate, observation times, etc.) and/or methodologies for evaluating motility inhibition of the worms [e.g., manual vs automated and quantitative motility measurements by the WormAssay (52) used in this study].The larval developmental assay has been used to measure drug efficacy because the failure to molt is likely linked to early establishment and development within the host and therefore clearance of the infections (53).The failure of larval development in H. contortus after staurosporine treatment demonstrates the broad-spectrum activity of this compound against the ruminant GIN parasites and suggests that additional in vitro screening assays should be performed to comprehensively evaluate drug efficacy.
In summary, we were able to validate the predicted broad-spectrum efficacy of K-252a, staurosporine, and digitoxigenin against ruminant GIN parasites in vitro and confirmed repositioning as a viable approach to developing broad-spectrum anthelmin tics, although improving treatment efficacy remains important.Combining these drugs into a single treatment may prove efficacious and may enable reductions in each drug's dosage, which deserves empirical confirmation.

DISCUSSION
Omics-driven drug target identification, combined with chemogenomic analyses, has previously identified drugs with anthelmintic potential proven through subsequent experimental confirmation (54)(55)(56)(57).Given the massive genetic diversity within the Nematoda, we chose to prioritize drug targets and compounds, using this approach, on a select taxonomic group: the GINs in the Trichostrongyloidea superfamily (clade Va), which commonly infect ruminants and pose the greatest economic challenge to producers worldwide (58)(59)(60)(61).To facilitate studies on the major species of veterinary importance, we undertook a stepwise approach (Fig. 1): we (i) sequenced, assembled, and annotated the genomes of three clade Va species with no or limited prior omics resources, enabling the creation of a database of five major GIN species of ruminants, (ii) employed a comparative multi-omics approach to identify conserved targets that may prove essential for parasite survival among all clade Va species, (iii) identified drugs and drug-like compounds capable of interacting with those candidate target molecules, and (iv) used the computational evidence to prioritize drugs for experimental screen ing against xL3 of clade Va, evaluating their repurposed function as broad-spectrum anthelmintics.
The GINs in the Trichostrongyloidea (clade Va) that were selected for genome sequencing, assembly, and annotation are globally important parasites of ruminants, and their polyparasitism (concurrent infections with multiple species) is common and leads to worse outcomes in livestock.Unfortunately, genomic resources for comparative studies and exploring aspects of parasite biology, evolution, and parasitism have been limited for these species.This has also curtailed the discovery of conserved targets for drug development.In this study, we newly assembled draft genomes for three ruminant GINs in clade Va (C.oncophora, T. colubriformis, and O. ostertagi).Based on the CEGMA assessment (when accounting for both complete and fragmented core eukaryote genes), our assembly completeness was 86.7%-91.5% and was within the range calculated for other clade Va species genomes (Table 1).A comparison of BUSCO completeness (v4, nematoda_odb10, n = 3,131) between the three new draft genomes and clade V species genomes available in WormBase Parasite indicated that our newly assembled genomes (76.1%-83.2%)closely align with the mean BUSCO scores (81.0%) estimated with the clade V genomes (Fig. S1).Furthermore, in the three new assemblies, over 90.0% of genes in each species were assigned to OGs identified in the seven clade V species, with more than 50.0% of the genes classified into VaC OGs (Fig. 3A).These results demonstrated that the quality of our new draft genomes is comparable to other previously published genomes of ruminant GINs in clade Va.Mitochondrial genome-based clustering has been performed previously (32), but this is the first whole genome-based clustering that inferred phylogenetic support for the subdivision of clade Va species (Fig. 2).The species tree inferred by OrthoFinder displayed two groups of ruminant parasites, clustered by taxonomy.By including outgroup host species (cattle, sheep, and goats), a root of the species tree was determined with a high STAG support value (99.9%), indicating a possible divergence of O. dentatum (Strongyloidea family) from a common ancestor of Trichostrongyloidea family (Fig. S7).Interestingly, branches for Trichostrongyloidea species had relatively low support values (between 10.1% and 23.3%) similar to those observed on a branch for small ruminant (58.9% for sheep and goats).We hypothesize that this is due to a recent speciation event among the target species resulting in insufficient phylogenetic information to resolve the order of terminal branches.The species tree inference may improve with more complete genome assemblies; however, these data generally support our overarching hypothesis that common drug targets exist for these closely related parasites.In addition, it has been suggested that predilec tion sites within hosts could drive modifications to the genomic sequence of parasites while adapting to habitat-specific challenges, which is then reflected in the divergence of evolutionary lineages (32,62,63).The species tree inferred with the newly assem bled genomes showed that among ruminant GINs in clade V, the overall clustering patterns confirmed the previously observed taxonomic classification (32-34): a possible speciation from a common ancestor of Trichostrongyloidea and Strongyloidea super families.Within the Trichostrongyloidea superfamily, the clade Va species residing in the abomasum/small intestine, where food digestion and absorption occur, clustered independent of the lungworm parasite, D. viviparus in clade Vb, suggesting the existence of a lineage-specific subset of genes related to the predilection sites within the hosts.
The taxonomic conservation and diversification of genes across the GIN of ruminants were investigated to identify drug target candidates presumably essential for the survival of clade Va parasites.Of great interest are the 220 VaSC OGs, where gene members selected based on lineage-specific conservation appeared to have minimal duplication events.Functional enrichment shows that genes among these 220 OGs have critical functions that may serve as promising drug targets.For example, the most enriched pathway for this group of OGs is "cellular senescence" (map04218; P ≤ 0.0013), which is related to the arrest of cellular development and proliferation (Fig. 3D).In most of the major trichostrongylid species, survival and development are influenced by environmen tal cues such as temperature and moisture.Early in the free-living stages, Teladorsagia and Trichostrongylus spp.larvae can long endure in the environment and develop throughout the year, even during the cold winter season, by slowing down their growth rate (64,65).In addition, most trichostrongylids undergo developmental arrest within the host in a season-dependent manner to increase pasture survival when environmental conditions become more permissive (59,66,67).Therefore, functional enrichment in cellular senescence may be related to the versatility of trichostrongylid species that modulate their developmental cycles for adapting to dynamic environments.Conse quently, those genes become interesting targets for anthelmintic drug discovery and development.
Comparisons of the transcriptional profiles of the genes in these 220 VaSC OGs identified three OGs with gene members overexpressed in parasitic compared to free-living stages (Fig. 5C).One of them, OG0009146, is predicted to be "α-crystal lin/hsp20." Functional roles of HSP20, especially in helminths, have been implicated in various cellular processes such as stress response, development, pathogenicity of parasites, host immune response, and nutrient uptake, where the worms must respond to sudden changes between environments and hosts and during adaptation to their habitats inside of the hosts (45,68,69).In this study, a dendrogram of all hsp20 genes identified in the seven clade V species and the non-parasitic nematode, C. elegans, determined two distinctive groups of "α-crystallin/hsp20" genes: one specific to the parasites and the other common to all nematode species examined including C. elegans (Fig. S4).In particular, the six hsp20 gene members, characterized by the highly conserved α-crystallin domain, in OG0009146 were located at one of the clade Va-spe cific branches, consistent with a descent from the last common ancestor of the intestinal parasites infecting the ruminant hosts (Fig. S4).Although further studies will be needed to characterize their functional roles, these observations raise the possibility that the HSP20 proteins in OG0009146 evolved to promote tolerance to cellular stress induced under hostile environments, thereby increasing the chances for survival within the host.It is worth noting that the RNA-seq data for C. oncophora and O. ostertagi, sourced from a prior publication, and the T. colubriformis data, generated in this study, each had a single replicate representing different life cycle stages (70).This potential limitation could impact the detection of DEGs when comparing various stages within each species.Nevertheless, the primary objective of our analysis was to identify potential drug target candidates conserved within the clade Va species.To achieve this, we adopted an approach that involved integrating gene-level expression patterns (CEGs/DEGs) in each species into OG-level expression patterns (CEOs/DEOs) that are conserved across the five clade Va species.This method was designed to minimize the likelihood of selecting genes that were not consistently expressed across the five clade Va species, allowing us to identify conserved transcriptional regulation for 220 VaSC OGs.In conclusion, the comparative genomic and stage-specific transcriptomic analyses revealed interspecies conservation of a subset of genes (proteins) and their potential relationships to the survival of clade Va parasites, implicating them as ideal targets for drug intervention.
Chemogenomics identified 11 compounds spanning targets across nine OGs (Table 2); some of them were known to have motility-inhibitory effects against nematode species, and others were newly proposed as anthelmintics in this study.One of the target classes identified was kinases (OG0009155, OG0009699, and OG0008885); the corresponding drug candidates that likely bind to those protein targets are stauro sporine, K-252a, calcitriol, and digitoxigenin.Both staurosporine (CHEMBL388978) and K-252a (CHEMBL281948, an analog of staurosporine) are non-selective, potent kinase inhibitors that suppress not only tumor cell growth but also parasitic infections (49)(50)(51)(71)(72)(73)(74)(75)(76).Staurosporine, a broad-spectrum anthelmintic, effectively inhibits the motility of phylogenetically diverse nematode species such as Trichuris muris (clade I), Brugia pahangi (clade III), Ascaris suum (clade III), and H. contortus (clade V) (49)(50)(51).K-252a interrupts host cell invasion of Trypanosoma cruzi by targeting TrkA kinase activity on host cells required for infection.It was also shown to directly interfere with calciumdependent protein kinase 1 activity in Plasmodium falciparum (74,75).Calcitriol, an active form of vitamin D as well as a steroid hormone, plays an important role in Ca 2+ homeostasis and binds to the membrane-associated vitamin D receptor (VDR), leading to rapid modulation of kinase signaling pathways such as those involving protein kinase C, MAP kinases, and calcium/calmodulin kinase II (77)(78)(79).In addition, calcitriol can elicit anti-tumor effects by inhibiting cellular proliferation and activating the apoptotic pathway by binding to the nuclear VDR (80)(81)(82).Although calcitriol is not known to be a kinase inhibitor, it was selected as a candidate lead compound based on drug-like properties against calcium/calmodulin-dependent kinase as well as its FDA approval status (Table 2).Digitoxigenin, one of the three cardiac glycosides along with digitoxin and digoxin, is an inhibitor of Na+/K + ATPase (sodium-potassium ATPase), a signal transducer converting extracellular signals into the activation of protein kinase cascades (83)(84)(85).This was selected based on potential drug-like efficacy on casein kinase II (pChEMBL 8.1 and QED 0.69; Table 2).
In vitro screening of the 11 prioritized drugs demonstrated, in general, that kin ase inhibitors are capable of causing significant motility inhibition or morphological abnormalities against GINs in clade Va.Considering the potency of staurosporine as a non-selective kinase inhibitor, its broad-spectrum activity against other nematode species was not surprising; however, here, we extend its application to affecting the motility of clade Va species.The IC 50 assay indicated that K-252a, the derivative of staurosporine, better inhibits the motility of GINs in clade Va.Most importantly, digitoxigenin was newly discovered as active against the clade Va species, but its efficacy was selective to GINs of the small intestines.Some of the compounds used for pathogenic infections such as degrasyn, anabaenopeptin B, and camptothecin did not show significant inhibitory effects on the motility of C. oncophora, T. colubriformis, or O. ostertagi; however, success in this assay was based primarily on phenotypic screenings of xL3 worms and was predicated upon movement changes or shedding of the old cuticles that are well-established and easily monitored (49)(50)(51)55).Other readouts may have proven more successful for these and other drugs.Therefore, the success rate of drug discovery may increase by examining other phenotypic changes (apoptotic cell death or tissue damage) or monitoring additional parasitic stages for changes in feeding, movement, development, and/or fertility of adult worms all of which could reduce morbidity and/or mortality caused by parasite infections.
Overall, these results demonstrate that a proof-of-principle study from knowledgebased identification of drug targets to chemogenomic prioritization of compounds inhibiting those targets is a powerful approach that engenders a success rate much higher than high-throughput phenotypic screening (86,87).In addition, the genomic and transcriptomic resources generated in this study could be used for improving diagnostics, therapeutics, and vaccines, along with population genomic studies that require new reference genomes to undertake genome-wide association studies to better understand drug resistance in these economically important parasites.

Parasites
To propagate parasites for DNA and RNA isolation or for producing exsheathed L3 (xL3) for in vitro motility inhibition assays, 1-day-old Holstein bull calves were transferred to indoor, concrete-floor pens to prevent exposure to parasites.Over a period of 8 weeks, all calves were provided water ad libitum and weaned from milk replacement to a diet of 16% protein supplemented with alfalfa and orchard grass hay.At 12 weeks of age, animals were pre-treated with fenbendazole (label dose) and then 4 days later, monospecifically infected with 75-200,000 infective L3 of C. oncophora (Ghent University), T. colubriformis (Zoetis Inc.), O. ostertagi (Beltsville MD strain), or H. contortus (strain UGA2004).Daily fecal collections were performed between 16 and 35 days post-inocu lation and cultured at room temperature in the presence of sterilized vermiculite and water for a period of 14 days.Infective L3 were collected using a Baermann funnel and secondarily purified over a fine mesh screen.Adult worms were collected by allowing infections to proceed for 21-25 days, then sacrificing the animals and collecting adult worms from either the abomasum (T.colubriformis, O. ostertagi, and H. contortus) or the upper 20% of the small intestines (C.oncophora) and purified over a Baermann funnel.Fourth-stage larvae (L4) were collected in a manner similar to adult worms except that infections were terminated at 9-10 days after infection.The purity of all preparations was evaluated by PCR as described (88,89).L3 of C. oncophora, T. colubriformis, and O. ostertagi were frozen at −80°C for later genomic DNA extraction, and L3, L4, and adult T. colubriformis were used immediately for RNA isolation; L3 of C. oncophora, T. colubriformis, O. ostertagi, and H. contortus were stored at 4°C until needed.To generate the xL3, infective L3 were incubated in 0.15% sodium hypochlorite at 37°C for 20 min, followed by five washes in 1× phosphate-buffered saline and then used immediately for in vitro motility inhibition assays with or without drug treatment.

Genome sequencing, assembly, and annotation
Genomic DNAs (gDNA) were extracted from 150,000 infective L3 treated with proteinase K:SDS at 65°C and supplemented with 100 mM β-mercaptoethanol, followed by organic extraction and ethanol precipitation.All samples were subsequently treated with RNase, then extracted and repurified as described above.For larger size DNA, the L3 were treated with Proteinase K:SDS:β-mercaptoethanol at 65°C for 2-4 min, then supplemen ted with molten 2% low melting point agarose, and the mixture was transferred to plug chambers and allowed to harden.The plugs were removed from the chambers and left to incubate in the digestion medium at 50°C overnight.The incubation was repeated twice with a fresh digestion medium.After extensive washing, the plugs were treated with agarase, and the resulting high molecular weight DNA was purified through multiple ethanol precipitations.

C. oncophora and O. ostertagi
Short insert (450 bp), 3 kbp, and 8 kbp whole genome shotgun libraries were gener ated and sequenced on the Illumina HiSeq 2000 sequencer, and PacBio subreads were generated on the PacBio RS machine.The Illumina whole-genome sequencingdata were assembled using AllPaths-LG, which reported a coverage depth of 91.7× for C. oncophora and 86.7× for O. ostertagi (90,91).The PacBio data were then used to fill gaps in the AllPaths-LG assembly using PBJelly (92).Nanocorr was used for error correction of the assembly (93).The Illumina paired-end reads were used to close gaps and extend contigs using the gap closure tool Pygap, which relies on the Pyramid assembler, a tool developed in-house at the McDonnell Genome Institute (MGI).The final polishing step was to run the L_RNA_scaffolder, which used RNA sequencing reads to improve contiguity (94).A repeat library was generated using Repeatmodeler (95), ribosomal RNA genes were identified using RNAmmer (96), and transfer RNAs were identified with tRNAscan-SE (97).Non-coding RNAs, such as microRNAs, were identified by a sequence homology search against the Rfam database (98).Repeats and predicted RNAs were then masked using RepeatMasker (95).Protein-coding genes were predicted using a combination of the ab initio programs Snap (99), Fgenesh (100), and Augustus (101), as well as the annotation pipeline tool Maker v2.31.10 (102), which aligns mRNA, expressed sequencing tags (ESTs), and protein information from the same species or cross-species to aid in gene structure determination and modification.A consensus gene set from the above prediction algorithms was generated using a logical, hierarchical approach developed at the MGI (103).Sequence Read Archive (SRA) accession numbers for all assembled genome sequence data and raw genomic reads are provided in Tables S1 and  S2.

T. colubriformis
10× genomics data were generated from L3-derived genomic DNA on the Illumina HiSeq X platform.The Supernova assembler v2.0.1 was run with parameters configured to down-sample the input data to approximately 50× coverage of the T. colubriformis genome (104).Initial assembly output was prepared using the supernova mkoutput function with a 1 kb contig size cutoff.PacBio Sequel subreads were generated as a supplement to the 10× Genomics data and used to scaffold the Supernova contigs using the SSPACE-LongRead program (105).We then ran the GapFiller program, which further improved the assembly by using the original 10× Genomics data to fill gaps in the scaffolded sequence (106).We ran multiple iterations of GapFiller until we stopped seeing improvement.PBJelly ( 92) "minmapq50 all q30" was run as well, and while it did not result in an improved assembly, it did capture 10 additional genes as detected by BUSCO (107).Those 10 genes, each buried within 5 kb flanking regions, were merged into the main assembly.Finally, we applied a polishing step using PILON (108).This tool uses the alignment of the original 10× Genomics reads back to the assembly to identify and correct sequencing errors.This produced our finished assembly.A repeat library was generated by running Repeatmodeler v1.0.11 on a 10 kb size filtered version of the finished assembly (95).This size filtration step resulted in a 70% reduction in the number of contigs while retaining 80% of the total assembly length and improved the quality of the final repeat library.We then generated a soft-masked version of the fully finished assembly using RepeatMasker (95).Braker2 (v2.1.0)(109) was then used to train Augustus gene models (110) for the downstream Maker execution (102).The soft-masked assembly and the collection of T. colubriformis RNA-seq listed in Table S3 were used for this training.Maker v2.31.10 was then run to generate the final gene set (102).Maker was given the fully finished assembly and allowed to do the repeat masking internally using the repeat library generated by Repeatmodeler (95) described above.Evidence provided to Maker included assembled transcripts from the previ ously mentioned RNA-seq samples generated using HiSat2 (v2.1.0)(111) and String Tie (v1.2.4) (112), the Augustus gene models built by Braker2, and protein sequence from nine related organisms retrieved from Wormbase Parasite (version WBPS14) (113): H. contortus (PRJEB506), T. circumcincta (PRJNA72569), O. dentatum (PRJNA72579), D. viviparus (PRJNA72587), Ancylostoma caninum (PRJNA72585), Ancylostoma ceylanicum (PRJNA72583), A. duodenale (PRJNA72581), Heterorhabditis bacteriophora (PRJNA13977), and Necator americanus (PRJNA72135).Maker was run using the setting keep_preds = 1 to enable the rescue of weakly supported genes in cases where InterProScan (114) detected a Pfam domain.Finally, gene product naming was performed using the PANNZER web tool (115).

Orthogroup inference and comparative genomic analyses
Nucleotide sequences were retrieved from the seven clade V species belonging to the Strongylida mentioned above.For all sequences except H. contortus, Redundans v0.14a was used to remove redundant contigs (possibly representing uncollapsed haplotypes) from the current version of assemblies without gap-closing (--nogapclosing) and scaffolding (--noscaffolding) steps (116), and the longest isoform of each gene was selected.The evaluation of annotated gene set completeness utilized BUSCO v4 (nematoda_odb10, n = 3,131) (107) and CEGMA v2.5 (n = 248) (117).For comparative purposes, 28 genomes of clade V species, obtained from WormBase Parasite, were downloaded and analyzed using BUSCO v4.The protein sequences were provided as an input to Orthofinder to infer orthogroups (OGs) in the five clade Va species with additional two outgroup species, D. viviparus (clade Vb) and O. dentatum (clade Vc), using the default settings (31).Functional annotation of genes in each species was predicted using InterProScan v5.28-67.0(114), GhostKOALA v2.2 (118), and BLAST+ v2.12.0 (119) search against MEROPS v12.1 (120).Putative functions of OGs were assigned if greater than or equal to 50% of gene members for each OG agreed on the corresponding functional annotation(s).Functional enrichment analyses of KEGG pathway and InterPro domain on subsets of OG such as VaS, VaC, and VaSC OGs were performed using over-representation analysis provided by WebGestalt (121) with a minimum of two OGs per term and P ≤ 0.05 threshold for significance.Computational Analysis of gene Family Evolution was used to model significant gene gain and loss in the seven clade V species (122).To investigate the gene-family dynamics, both an ultrametric tree generated from a rooted species tree (as produced by OrthoFinder) and a subset of OGs, where the number of gene members was less than 100, were taken into consideration.The functional enrichment analysis on the rapidly evolving OGs was performed using WebGestalt with the same criteria mentioned above.The number of rapidly evolving proteases and their relative ratios to aspartic, cysteine, metallo-, or serine proteases identified by MEROPS within each species was visualized along with the phylogenetic tree using Evolview (123).Protein sequence data of three host species (Bos taurus: GCF_002263795.2, Ovis aries: GCF_016772045.1, and Capra hircus: GCF_001704415.2) retrieved from NCBI Genome (124) were included in additional OrthoFinder analyses along with the seven clade V parasites to identify parasite genes orthologous to host protein-coding genes.

T. colubriformis RNA-seq production
L3, L4, and adult T. colubriformis were used fresh.Approximately 75-100 μL of settled worms were disrupted in a Dounce homogenizer in the presence of 1.5 mL of Trizol.The resulting RNA was DNase treated and column purified (Zymo Research) to remove contaminating gDNA.cDNA libraries were prepared from RNA samples with random primers, and processed cDNA was sequenced on the Illumina NovaSeq 6000 platform (paired-end 150 bp reads).

Comparative transcriptomic analysis
RNA sequencing data used in this study include T. colubriformis libraries generated above and C. oncophora and O. ostertagi libraries deposited in the Sequence Read Archive NCBI database (70) (Table S3).RNA sequencing reads were trimmed for adapters using Trimmomatic v0.36 (125), mapped to their respective genomes either listed in Table S1 or mentioned above using STAR v2.6.1d(126), and quantified per gene using featureCounts v1.5.2 (127).The read counts were combined for technical replicates when applicable.
Based on the availability of RNA-seq data within each species and stage, the life cycles of the parasites were categorized into three groups/stages: free-living (from egg to sheathed L3), parasitic (xL3, L4, and adult), and all stages.Comparative tran scriptomic analyses were performed in two ways at the gene level followed by an OG level: three-groups comparison (free living, parasitic, and all) for constitutive expression analysis and two-groups comparison (free living and parasitic) for differential expression analysis.Transcriptional conservation of OGs across species was examined for VaSC expression (220 VaSC OGs), where at least one gene member from each of the clade Va species exhibited expression in one of the stages defined above.
Constitutively expressed genes within each species were identified using normalized gene expression values (FPKM) calculated for each gene in each stage."Constitutive" expression of the gene across all free-living and/or parasitic stages was defined as requiring FPKM ≥ 1 at all available corresponding stages of the parasite.The 1 FPKM cutoff was selected based on it being used in previous literature (128)(129)(130) and an initial study that identified it as a cutoff at which "transcript detection was robust" for a typical 2-kb mRNA (131).Constitutive expression orthogroups were then identified based on the membership of CEGs spanning all five clade Va species within each OG and were classified as either CEOs in "free living, " "parasitic, " and/or "all." Differential gene expression analysis was performed using transcriptomic evidence within each species.Significant upregulation of genes in either free-living or parasitic stages was evaluated using DEseq2 v1.34.0 with absolute log 2 fold change ≥ 1 and a threshold for significance P ≤ 0.05 (132).Of the 220 VaSC OGs, differential expression orthogroups were identified with conserved gene expression patterns that were higher in either free-living or parasitic stages across all five species.
OrthoFinder (31) was also run separately to identify the heat-shock protein 20 ortholog of the gastrointestinal parasites in clade V, based on orthology to the non-para sitic nematode, C. elegans (PRJNA13758.WS277), obtained from Wormbase Parasite (113).The protein sequences of all hsp20 orthologs identified were used to perform multiple sequence alignment by ClustalW (133) with BLOSUM weight matrix in a default setting followed by converting the alignment file to a nexus format using ALTER (134).The maximum likelihood phylogenetic tree of the hsp20 genes was then constructed by RaxML v8.2.12 (135) using "PROTGAMMAWAG" substitutional model with a bootstrap value of 1,000 and visualized by "ggtree" R package (136).

Chemogenomic screening using the ChEMBL database
To identify drug-like compounds that can be repurposed as anthelmintics against the ruminant GIN in clade Va, a BLASTp (119), search with e-value ≤ 10 −4 cutoff was per formed based on sequence similarity between 2,221 target proteins in the 220 VaSC OGs and target proteins deposited in ChEMBL v27 (47).Refinement of the druggable targets was performed based on the overlap of the significant target hits found in all five clade Va species for each OG.The identified target hits were linked to ChEMBL compounds, where the corresponding binding affinity data are available from the drug:target affinity score (pChEMBL ≥ 5), which resulted in one-to-one or one-to-many drug-target pairs per OG.Expression profiles of 220 VaSC OGs assessed in the previous step were also considered when selecting parasite target gene products transcribed in all clade Va species.The ChEMBL drug-target pair(s) of each OG, if applicable, were then ranked based on weighted QED, pChEMBL scores, commercial availability provided by the ZINC database (48), and cost (Tables S7 and S8).

In vitro screening of motility of parasites
Approximately 80 freshly generated xL3 were placed in each well of a 96-well plate and cultured in 50 µL of the Luria-Bertani medium containing 100 units/mL of penicillin, 100 µg/mL of streptomycin, and 250 µg/mL of amphotericin B at 37°C with 5% CO 2 for C. oncophora and T. colubriformis and 20% CO 2 for O. ostertagi and H. contortus.A list of drugs used in the in vitro assays and relevant information (sources, solvents, controls, stock, and test concentrations) are summarized in Table S8.For each drug treatment, three biological replicates were conducted.
The motility of the worms was recorded for 20 seconds daily up to 3 days using a Keyence BZ-X810 microscope with 2× magnification.WormAssay v1.7.1 (52) with the "consensus voting luminance difference" algorithm was used to quantify the movement of the worms in each replicate and enumerate using an arbitrary unit (mean movement units, MMUs) (52).Relative motility inhibition (%) was calculated by dividing the MMUs of the treated worms by the average MMUs of control worms treated with the corre sponding solvent such as 1% DMSO or methanol, subtracting this value from 1, and then multiplying by 100 (52).Drug efficacy was evaluated according to the following criteria: ≥75% motility inhibition relative to the controls as "effective" and ≥50% motility inhibition as "moderate." After the initial in vitro screening at a single concentration, IC 50 assays were performed using six-point dilutions (from 100 to 0.3 µM) with three biological replicates for the following active drugs: digitoxigenin, K-252a, and staurospor ine.The IC 50 values of the compounds were calculated using Prism v9 with a non-linear regression curve (40,52).For in vitro treatment of the three active compounds on H. contortus, xL3 were incubated at 37°C with 20% CO 2 , and the presence of empty cuticles was observed to estimate their development to L4 at days 0, 1, 2, 3, and 7 using a Keyence BZ-X810 microscope with the 20× magnification.

DIRECT CONTRIBUTION
This article is a direct contribution from Makedonka Mitreva, a Fellow of the Ameri can Academy of Microbiology, who arranged for and secured reviews by Isheng Tsai, Academia Sinica, and James Cotton, University of Glasgow.

FIG 1
FIG 1 Overview of omics-driven drug discovery workflow used in this study.(A) Genomes of key clade Va species of interest (indicated by green arrows) were newly assembled and annotated, providing a critical component of the comparative analysis.(B) Drug target prioritization was performed based on the identification of orthologous genes (orthogroups) and selection of 220 VaSC OGs that were exclusively conserved in the interest group of species (clade Va).Stage-specific transcriptomic profiles were examined to functionally characterize drug target candidates spanning the species of interest.(C) Compounds targeting prioritized drug candidates from panel B were predicted based on comparing their protein sequence similarity with ChEMBL drug targets, followed by exploiting binding affinity to the ChEMBL drug targets, commercial availability, etc. (D) In vitro screening of exsheathed L3 (xL3) were used to validate their predicted efficacy on the species of interest.Figures were created with Biorender.com.

FIG 3
FIG 3 Identification, classification, and functional enrichment of orthogroups in clade V species.(A) A total of 17,689 OGs identified in clade V species were classified into three groups: clade Va-conserved (VaC) OGs, clade Va-specific (VaS) OGs, and clade Va-specific and -conserved (VaSC) OGs.The percentage (%) of assigned genes to each of the OG classifications was calculated per species.In addition, the percentage of genes in OGs shared by D. viviparus or O. dentatum (non-clade Va species), referred to as non-clade Va-specific OGs, was evaluated for the five clade Va species.The percentage of genes in species-specific OGs, consisting entirely of genes from one species, was assessed per species.(B-D) Significant KEGG pathway enrichment in 6,486 VaC OGs (B), 4,706 VaS OGs (C), and 220 VaSC OGs (D).

FIG 4
FIG 4 Overview of comparative transcriptomic analysis across the five ruminant GIN in clade Va.Transcriptional levels of all genes in each species were evaluated at a gene level with FPKM values and DESeq2 results.For an OG level analysis, conserved expression patterns of gene members in the 220 VaSC OGs were assessed across the five clade Va species.Accordingly, the 220 VaSC OGs were assorted to one of the following categories: CEO in all, CEO in free-living, CEO in parasitic, DEO in free-living, and DEO in parasitic stages.

FIG 5
FIG 5 Identification of stage-specific transcriptomic profiles and conserved expression patterns among ruminant GINs in clade Va.(A) Within each species, constitutive expression of genes in all, free-living, and parasitic stages were examined, and their relative frequencies were represented in a bar graph.(B) Same as panel A but for differentially expressed genes: significantly higher in free-living and parasitic stages.(C) Conserved expression patterns at the OG level were determined based on the presence of at least one representative gene being expressed in each of the clade Va species and agreement on one of the five expression categories in all clade Va species.The pie chart shows a summary of the expression patterns of 220 VaSC OGs.

FIG 6
FIG6 Overall workflow to identify and prioritize drug-like compounds having anthelmintic potentials to the ruminant GINs (clade Va species).Homology relationship was inferred between 2,221 parasite genes (220 VaSC OGs) and 9,211 drug targets in ChEMBL using Blast sequence similarity searches.The OG-level integration of the significant blast hits resulted in 882 consensus targets associated with 40 OGs.Next, ChEMBL target:compound pair information such as pChEMBL and expression profiles of the 220 VaSC OGs were considered to determine 57,460 compounds that could be repositioned as anthelmintics.After selecting up to 20 drug-like compounds per OG, followed by assessing commercial availability using the ZINC15 database, a total of 11 compounds potentially associated with nine OGs were identified, and their efficacies were evaluated in vitro.

FIG 7
FIG 7 In vitro screening of three kinase inhibitors on exsheathed L3 (xL3) of gastrointestinal parasites in ruminants.Movement of xL3 C. oncophora, T. colubriformis, and O. ostertagi was observed after in vitro treatment of digitoxigenin (A), K-252a (B), and staurosporine (C) for 3 days.Control worms were treated with either 1% DMSO or 1% methanol according to a solvent used for each drug dilution.The percentage of motility inhibition was calculated relative to control worms.For each treatment, three biological replicates were prepared.(D) IC 50 of digitoxigenin, K-252a, and staurosporine on day 3.

TABLE 1
Genome assembly and annotation statistics of C. oncophora, T. colubriformis, and O. ostertagi with other major ruminant parasites in clade V (the order Strongylida)

C. oncophora T. colubriformis O. ostertagi H. contortus T. circumcincta D. viviparus O. dentatum
Gene statistics based on the presence of start and stop codons and the absence of internal stop codons.
a b Annotated genes inferred by Gene Ontology, InterPro domain, and KEGG pathway.

TABLE 2 A
list of 11 prioritized compounds for nine VaSC OGs a,b a Constitutively expressed orthogroup refers to orthogroup having a membership of constitutively expressed genes spanning all the five clade Va species in all, free-living, or parasitic stages.b Differentially expressed orthogroup refers to orthogroup having a membership of differentially expressed genes spanning all the five clade Va species between free-living and parasitic stages.