Introduction
Haemonchus contortus (barber’s pole worm), is one of the most economically important pathogenic nematodes that attacks small ruminants, such as sheep and goats, representing a global animal health issue through drastic losses in livestock. The life cycle of this nematode is comprised of six stages, namely, eggs, first to fourth stage larvae (L1, L2, L3 and L4) and adults (Veglia, 1915). In Haemonchus contortus and related nematodes, infective L3s are ingested by the host and enter the rumen where larval exsheathment occurs that marks the transition from the free-living to the parasitic stages of these parasites (Rogers, 1960). Molecular and omics driven studies targeting this key transition phase can help us understand the developmental physiology of this model nematode which can be used to develop both biological and biotechnological control strategies in the future.
There have been earlier transcriptomic and expressed sequence tag (EST) studies describing the differences in transcription between free-living infective (iL3) and parasitic (xL3) third-stage larvae of H. contortus (Cantacessi et al., 2010; Laing et al., 2013; Schwarz et al., 2013). However, a major knowledge gap remains amongst the available RNA-seq data sets for this parasitic nematode (Jex et al., 2019). However, all these studies applied the common laboratory practice of the utilization of sodium hypochlorite as the desheathment agent. Recently, a closed in vitro parasite culture system that effectively mimics rumen conditions to effectively stimulate exsheathment without chemical interventions has been reported (Bekelaar et al., 2018; Palevich et al., 2022; Palevich et al., 2023a). Briefly, this system involves an increase to rumen temperature (39°C) and a strictly anaerobic environment of predominantly carbon dioxide (CO2). The larval exsheathment method used in this study for an RNA-seq approach, has recently been validated as an adaptable technique for multi-omic’ applications of H. contortus (Palevich et al., 2022). The aim of our study was to provide the scientific community a detailed time-series transcriptomic description of natural larval exsheathment of H. contortus, a valuable resource for future multiomics functional studies investigating larval exsheathment and parasite development.
Value of the data
• Time-series RNA-seq dataset using an in vitro system that effectively reproduces the natural rumen environment fills the remaining gap of the H. contortus developmental transcriptome.
• Three distinct expression profiles revealed post-trigger application and identification of the exact time interval that marks the rapid and irreversible population shift in larval exsheathment.
• The presented transcriptomic resource can guide future gene discovery research and the identification of novel compounds with broad-spectrum efficacy that can be used to either trigger premature exsheathment or inhibit the process all together.
Data
Haemonchus contortus is a serious threat in small ruminants causing serious animal production issues worldwide. To improve our understanding of the fundamental genetics of larval exsheathment, the transcriptomes of 100 samples of H. contortus L3’s at different stages of the in vitro exsheathment process (i.e., 49 time points over 24 h) have been sequenced (RNA-seq) using the Illumina NovaSeq6000 technology. We applied an RNA-seq approach to conduct a time-series transcriptome profiling of infective larvae pre- and post-exsheathment trigger application of the barber’s pole worm. Samples were collected, RNA isolated using TRizol method and stranded mRNA sequencing was carried out using the Illumina NovaSeq6000 platform (Macrogen, Inc.). RNA integrity number (RIN) for all samples were >7 with a total number of 5,918,887,284 clean reads generated. Duplicate runs were performed for all time points (except pre-treatment and t = 0 samples that were performed in triplicate) and raw sequences were submitted to Sequence Read Archive (SRA) with data deposited under the GenBank BioProject accession number PRJNA517503. Detailed information regarding BioSample and SRA accession numbers and other related statistics are provided in Supplementary Table S1. Transcripts were mapped against the reference H. contortus NZ_Hco_NP v1.0 genome (Palevich et al., 2019b; Palevich et al., 2023b) using STAR v.2.7.1a (Dobin et al., 2013).
The transcriptome data has revealed significant differences in the overall gene expression profiles of the parasite populations at different stages of the exsheathment process. A comparison of the total numbers of sequencing read pairs at each time point reveals a dramatic decrease at the 28-min interval (Figure 1B). This transition indicates a strong negative correlation between the total numbers of exsheathed H. contortus L3 larvae (xL3) and sequencing read pairs observed with respect to the total population of exsheathed xL3’s over time for these samples, that was strongest at the 28-min interval indicating the rapid and irreversible population shift in exsheathment.
FIGURE 1. Overview of the study design and experimental procedure. A closed in vitro system that effectively reproduces the two basic components of an anaerobic rumen environment (CO2 and 39°C) was used to trigger exsheathment (xL3) in Haemonchus contortus third-stage infective larvae (iL3) in O2-free CO2 saturated saline solution (A). Scanning electron micrographs of the exsheathment process with scale bars adjusted according to magnification of each image. Time-series RNA-seq analysis of Haemonchus contortus larval exsheathment (B). Comparison of the total numbers of clean sequencing reads versus total population exsheathment of L3 larvae. Total reads represent the sum of both unique and duplicate read pairs. Mean of the total exsheathment percentage (±SEM) at each time point across replicates. Top 100 differentially transcribed genes (DTGs) post exsheathment trigger application (C). Standardised relative abundances of inferred gene transcripts (Value) indicating low to high (red to green) relative gene expression (logFC≥2, FDR≥0.05) at each time point, and single linkage clustering. Three distinct profiles of expression were observed within the time-series RNA-Seq data: ‘prolonged/constant’ (pink), ‘oscillating’ (cyan, ≤25 min) and ‘ramped’ (purple, ≥28 min) expression for the majority of time intervals.
Overall, a total of 863 significantly (logFC≥2, FDR≥0.05) differentially transcribed genes (DTGs) post trigger application with three distinct patterns or profiles of expression observed (Figure 1C). A large proportion of these DTGs are predicted to encode hypothetical proteins or proteins that are yet to be characterized (Palevich et al., 2018), and genes that were only upregulated at the final 720 and 1,440 min (12 and 24 h) time intervals. Further, a significant decrease in the total sequencing read counts assigned to Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways was also observed at the 28-min interval (Figure 2). Across all time points, the categories ‘signal transduction’ and ‘carbohydrate metabolism’ were the two most abundant B-level KEGG functions (Figure 2B). The presented time-series RNA-seq dataset of natural larval exsheathment in H. contortus can be used for applications such as gene discovery or comparison with genomes and/or transcriptomes of other parasitic nematodes of veterinary importance to understand the molecular mechanisms driving their parasitism.
FIGURE 2. Time-series functional profiling of larval exsheathment in Haemonchus contortus. For each transcriptomic dataset, only the top level-A (A) and corresponding level-B (B) KEGG pathways are depicted. The total sequencing reads assigned to KEGG categories are averages across each time point.
Materials and methods
Nematode production and procurement
Pure cultures of H. contortus third-stage larvae (L3) were maintained by regular passage through five otherwise parasite-free lambs housed indoors at the AgResearch’s Grasslands campus (Figure 1A). L3 were cultured in fresh faecal material containing eggs collected into faecal bags on infected sheep. Faeces were pooled and mixed with vermiculite then placed in trays, moistened with tap water (at 20°C), covered and cultured for 10 days at 22°C–24°C. A modified Baermann technique was used to clean and separate the larvae from faeces (Viglierchio and Schmitt, 1983). Briefly, approximately 150 g of faeces was enclosed in paper facial tissues and suspended over a large conical measuring flask filled with unchlorinated water. The samples were incubated for 20 h and larvae washed by their movement through the apparatus. The faeces were then removed and the liquid carefully siphoned to a remaining volume of about 20 mL. The resulting L3s and volume were examined in a counting chamber (Whitlock S.F.E.L.O., Australia; volume 2 mL) at ×100 magnification and stored in tap water at 10°C until required. A single batch of larvae were used for all experiments to account for any batch-related confounding factors with larval viability and motility checked microscopically prior to any further experimentation.
The provenance of genomic DNA was verified with a 100% similarity identity to the representative chromosome-level genome of the anthelmintic-susceptible H. contortus field strain (Palevich et al., 2019a; Palevich et al., 2019b; Palevich et al., 2020), by automated Sanger sequencing of the second internal transcribed spacer (ITS-2) of nuclear ribosomal DNA following PCR amplification from genomic DNA.
In vitro larval exsheathment, sample collection and electron microscopy
Prior to in vitro testing, pooled L3 cultures were cleaned using autoclaved phosphate buffered saline (1× PBS) solution (137 mM NaCl, 2.7 mM KCl, 8 mM Na2HPO4, and 2 mM KH2PO4, pH 7.4) and acclimated overnight to room temperature by gravity migration filtration through nylon mesh (pore size 20 μm). Larval viability and motility were checked microscopically and quantified using a Petroff-Hausser chamber (Hausser Scientific) according to the manufacturer’s instructions.
In this study, we triggered larval exsheathment using a modification of the method described by Bekelaar et al., 2018, in a ‘closed’ in vitro system that effectively and reproducibly simulates the physiological conditions of the rumen (39°C and high CO2 concentration) as previously described (Palevich et al., 2022; Palevich et al., 2023a). Briefly, anaerobic 1× PBS solution was mixed in boiling dH2O and cooled to room temperature under a continuous flow of O2-free CO2. Once cooled, the PBS solution was transferred to 100 mL serum bottles in 70 mL aliquots and flushed with CO2 for 1 h. The bottles were sealed with butyl rubber bungs and aluminium crimp caps before being autoclaved at 121°C for 20 min. Larval cultures were transferred anaerobically (approximately 300,000 L3/sample) using a O2-free CO2-flushed 3 mL syringe and wide bore hypodermic needle (16 G thickness and length of 1–1/2 inch, BD), into aluminium-wrapped and pre-warmed serum bottles containing 70 mL of autoclaved PBS and incubated anaerobically at 39°C with gentle horizontal shaking at 75 rpm for up to 24 h in darkness.
Exsheathment of L3 larvae (xL3) was determined retrospectively by either complete or partial loss of the sheath and was measured by additional sub-sampling of each replicate beginning at t = 0 min up to 24 h after incubation at 39°C and CO2 anaerobic conditions. At each time point, contents were thoroughly mixed before 1 mL was transferred anaerobically to a 24-well plate and exsheathment enumerated. The numbers of xL3s in each subsample was quantified via 300-fold dilution of sample to another 24-well plate containing dH2O to yield approximately 300 larvae per replicate. Larvae were killed by the addition of 1 drop of 3% helminthological iodine solution (Lugol’s solution) and exsheathment enumerated. Samples were collected anaerobically from each time point, snap-frozen in liquid nitrogen and stored at −80°C until further processing.
Scanning electron microscopy (SEM) was performed as previously described (Palevich et al., 2019b; Palevich et al., 2022). Briefly, cryopreserved worms were gently spun, washed 3x in PBS and fixed in SEM primary fixative (3% glutaraldehyde, 2% formaldehyde in 0.1 M Phosphate Buffer pH 7.2) for 2 days at room temperature. Samples were dehydrated in a graded ethanol series, i.e. 25%, 50%, 75%, and 95% for 10–15 min each and 2× in 100% ethanol for 1 h, then Critical Point (CP) dried using liquid CO2 and mounted onto an aluminium specimen support stub using double-sided adhesive tape. Samples were sputter coated with gold (200 s) and observed using a FEI Quanta 200 Environmental Scanning electron microscope with energy dispersive x-ray spectroscopy (EDAX) module.
Total RNA extraction, library preparation and RNA-Seq
Snap-frozen samples containing 100 µL of packed worms were lysed using an 18 V drill loaded with a disposable RNAase-free polypropylene Micro-Pestle (Qiagen) until the mix is ground to a fine white powder. To the ground sample 250 µL of pre-warmed (40°C) TRizol (Life Technologies) was added and mixed thoroughly according to the manufacturer’s instructions, snap-frozen in liquid N2 and the homogenization of snap frozen samples in TRizol was repeated for five rounds in total to ensure complete disruption of the sample. To the homogenized sample, 750 µL of pre-warmed TRizol and 0.1 volume of chloroform were added, thoroughly mixed and centrifuged at 20,000 × g for 10 min at 4°C. The upper aqueous phase was transferred into a new Eppendorf tube and an equal volume of isopropanol and 0.1 volume of 3 M sodium acetate (pH 5.5) were added and gently mixed, and the mixture was stored at −20°C overnight. The RNA pellets were precipitated with ethanol, re-suspended in nuclease-free water (Life Technologies) and DNase I treated. RNA yield and quality were assessed using the Bioanalyzer 2,100 with the RNA 6000 Nano assay reagent kit from Agilent (Santa Clara, CA) and stored at −80°C. NGS sequencing libraries were generated from 1 μg of total RNA using TruSeq Stranded RNA Sample Prep Kit (Illumina) according to the manufacturer’s protocol. The resulting cDNA libraries were then paired-end sequenced (2 × 100 bp) using the NovaSeq6000 instrument (Macrogen, Inc.), to produce 2,959,443,642 read pairs in total with an average of 29,594,436 per sample.
Functional annotation and differential expression analysis of time-series RNA-seq data
Complete paired-end sequences were obtained and converted into FASTQ files (forward and reverse) using bcl2fastq (Illumina) package from Macrogen. Adaptor sequences, minimum length (36 bp), bad fraction (>10%) and low-quality bases with PHRED scores (Q) ≤ 20, were removed and paired-end reads cleaned using Trimmomatic v0.4 (Bolger et al., 2014). The reads were mapped against the reference of H. contortus NZ_Hco_NP v1.0 genome (Palevich et al., 2019b; Palevich et al., 2023a) using STAR version v2.7.1a (Dobin et al., 2013). The counts were compiled and tabulated, and a differential expression was performed using the DESeq2 v1.30.2 package (Love et al., 2014) with default parameters. A differentially transcribed gene (DTG) was defined as a transcript with an absolute log fold change ≥2 and an FDR of less than 0.05. The “top” DGTs were chosen by their ranking in absolute fold change and displayed a consistent level expression for the majority of time intervals. The sample features, GenBank BioSample and Sequence Read Archive (SRA) accession numbers for the time-series exsheathment transcriptomes of H. contortus are listed in Supplementary Table S1.
Possible protein coding regions within the genome-based transcripts were identified using the TransDecoder program v5.5.0 implemented in TRINITY v2.14.0 (Grabherr et al., 2011). The protein-coding regions were searched against the NCBI NR protein sequence database using the blastp function of DIAMOND v2.0.6 (Buchfink et al., 2015) with the output format of XML being specified. The results were imported into OmicsBox v1.4.11 (https://www.biobam.com), where the Blast2GO (Conesa et al., 2005) annotation functions were used with default settings. InterProScan v5.50–84.0 (Jones et al., 2014) and EggNOG-Mapper v1.0.3 with EggNOG v5.0.0 (Huerta-Cepas et al., 2018) were further used with default settings to annotate the predicted proteins. For analysis of KEGG biological pathways, gene abundance tables were generated by alignment of the sequencing reads to the NCBI NR database using the blastx function of DIAMOND v2.0.6 (Buchfink et al., 2015) aligner. The results were “MEGANised” using the tools provided with MEGAN6 Ultimate Edition (Huson et al., 2016) and loaded into the MEGAN software to assign putative functions to the DIAMOND alignment files for level A and B KEGG categories (Supplementary Table S2).
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.
Ethics statement
The animal studies were approved by the AgResearch Grasslands Animal Ethics Committee under the Animal Welfare Act 1999 in New Zealand (AEC application number 13928). The studies were conducted in accordance with the local legislation and institutional requirements. Written informed consent was obtained from the owners for the participation of their animals in this study.
Author contributions
NP: Conceptualization, Data curation, Formal Analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing–original draft, Writing–review and editing. PM: Formal Analysis, Methodology, Software, Validation, Writing–review and editing. RS: Funding acquisition, Supervision, Writing–review and editing.
Funding
The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This work was supported by the Agricultural and Marketing Research and Development Trust (AGMARDT) Postdoctoral Fellowship Programme (Grant P17001) awarded to NP.
Acknowledgments
We wish to acknowledge the helpful staff at the Manawatū Microscopy and Imaging Centre at Massey University, Palmerston North, New Zealand for assistance with conducting the electron microscopy.
Conflict of interest
Authors NP, PM and RS were employed by AgResearch Limited.
The reviewer IM declared a past co-authorship with the author NP to the handling editor.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2023.1257200/full#supplementary-material
References
Bekelaar, K., Waghorn, T., Tavendale, M., McKenzie, C., and Leathwick, D. (2018). Heat shock, but not temperature, is a biological trigger for the exsheathment of third-stage larvae of Haemonchus contortus. Parasitol. Res. 117, 2395–2402. doi:10.1007/s00436-018-5927-2
Bolger, A. M., Lohse, M., and Usadel, B. (2014). Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120. doi:10.1093/bioinformatics/btu170
Buchfink, B., Xie, C., and Huson, D. H. (2015). Fast and sensitive protein alignment using DIAMOND. Nat. Methods 12, 59–60. doi:10.1038/nmeth.3176
Cantacessi, C., Campbell, B. E., Young, N. D., Jex, A. R., Hall, R. S., Presidente, P. J. A., et al. (2010). Differences in transcription between free-living and CO2-activated third-stage larvae of Haemonchus contortus. BMC Genomics 11, 266. doi:10.1186/1471-2164-11-266
Conesa, A., Götz, S., García-Gómez, J. M., Terol, J., Talón, M., and Robles, M. (2005). Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics 21, 3674–3676. doi:10.1093/bioinformatics/bti610
Dobin, A., Davis, C. A., Schlesinger, F., Drenkow, J., Zaleski, C., Jha, S., et al. (2013). STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, 15–21. doi:10.1093/bioinformatics/bts635
Fu, L., Niu, B., Zhu, Z., Wu, S., and Li, W. (2012). CD-HIT: accelerated for clustering the next-generation sequencing data. Bioinformatics 28, 3150–3152. doi:10.1093/bioinformatics/bts565
Grabherr, M. G., Haas, B. J., Yassour, M., Levin, J. Z., Thompson, D. A., Amit, I., et al. (2011). Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat. Biotechnol. 29, 644–652. doi:10.1038/nbt.1883
Huerta-Cepas, J., Szklarczyk, D., Heller, D., Hernández-Plaza, A., Forslund, S. K., Cook, H., et al. (2018). eggNOG 5.0: a hierarchical, functionally and phylogenetically annotated orthology resource based on 5090 organisms and 2502 viruses. Nucleic Acids Res. 47, D309–D314. doi:10.1093/nar/gky1085
Huson, D. H., Beier, S., Flade, I., Górska, A., El-Hadidi, M., Mitra, S., et al. (2016). MEGAN community edition - interactive exploration and analysis of large-scale microbiome sequencing data. PLOS Comput. Biol. 12 (6), e1004957. doi:10.1371/journal.pcbi.1004957
Jex, A. R., Gasser, R. B., and Schwarz, E. M. (2019). Transcriptomic resources for parasitic nematodes of veterinary importance. Trends Parasitol. 35, 72–84. doi:10.1016/j.pt.2018.09.010
Jones, P., Binns, D., Chang, H.-Y., Fraser, M., Li, W., Mcanulla, C., et al. (2014). InterProScan 5: genome-scale protein function classification. Bioinformatics 30, 1236–1240. doi:10.1093/bioinformatics/btu031
Laing, R., Kikuchi, T., Martinelli, A., Tsai, I. J., Beech, R. N., Redman, E., et al. (2013). The genome and transcriptome of Haemonchus contortus, a key model parasite for drug and vaccine discovery. Genome Biol. 14, R88–R16. doi:10.1186/gb-2013-14-8-r88
Love, M. I., Huber, W., and Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550. doi:10.1186/s13059-014-0550-8
Palevich, N., Britton, C., Kamenetzky, L., Mitreva, M., De Moraes Mourão, M., Bennuru, S., et al. (2018). Tackling hypotheticals in helminth genomes. Trends Parasitol. 34, 179–183. doi:10.1016/j.pt.2017.11.007
Palevich, N., Maclean, P., Baten, A., Scott, R., and Leathwick, D. M. (2019a). The complete mitochondrial genome of the New Zealand parasitic roundworm Haemonchus contortus (Trichostrongyloidea: haemonchidae) field strain NZ_Hco_NP. Mitochondrial DNA Part B 4, 2208–2210. doi:10.1080/23802359.2019.1624634
Palevich, N., Maclean, P. H., Baten, A., Scott, R. W., and Leathwick, D. M. (2019b). The genome sequence of the anthelmintic-susceptible New Zealand Haemonchus contortus. Genome Biol. Evol. 11, 1965–1970. doi:10.1093/gbe/evz141
Palevich, N., Maclean, P. H., Choi, Y.-J., and Mitreva, M. (2020). Characterization of the complete mitochondrial genomes of two sibling species of parasitic roundworms, Haemonchus contortus and Teladorsagia circumcincta. Front. Genet. 11, 573395. doi:10.3389/fgene.2020.573395
Palevich, N., Maclean, P. H., Candy, P. M., Taylor, W., Mladineo, I., and Cao, M. (2022). Untargeted multimodal metabolomics investigation of the Haemonchus contortus exsheathment secretome. Cells 11 (16), 2525. doi:10.3390/cells11162525
Palevich, N., Maclean, P. H., and Cao, M. (2023a). Non-targeted multimodal metabolomics data from ovine rumen fluid fractions. Microbiol. Resour. Announc. 12 (7), e0039223–23. doi:10.1128/MRA.00392-23
Palevich, N., Maclean, P. H., Carbone, V., Jauregui, R., and Umair, S. (2023b). Multi-omic profiling, structural characterization and potent inhibitor screening of evasion-related proteins of a parasitic nematode, Haemonchus contortus, surviving vaccine treatment. Biomedicines 11 (2), 411. doi:10.3390/biomedicines11020411
Rogers, W. (1960). The physiology of infective processes of nematode parasites; the stimulus from the animal host. Proc. R. Soc. Lond. Ser. B. Biol. Sci. 152, 367–386. doi:10.1098/rspb.1960.0045
Schwarz, E. M., Korhonen, P. K., Campbell, B. E., Young, N. D., Jex, A. R., Jabbar, A., et al. (2013). The genome and developmental transcriptome of the strongylid nematode Haemonchus contortus. Genome Biol. 14, R89–R18. doi:10.1186/gb-2013-14-8-r89
Veglia, F. (1915). The anatomy and life-history of Haemonchus contortus. Reports of the Director of Veterinary Research.
Keywords: Haemonchus contortus, parasitic nematode, transcriptome, exsheathment, rumen
Citation: Palevich N, Maclean PH and Scott RW (2023) Time-series transcriptomic profiling of larval exsheathment in a model parasitic nematode of veterinary importance. Front. Cell Dev. Biol. 11:1257200. doi: 10.3389/fcell.2023.1257200
Received: 12 July 2023; Accepted: 31 October 2023;
Published: 13 November 2023.
Edited by:
Juan Jesus Tena, Spanish National Research Council (CSIC), SpainReviewed by:
Ivona Mladineo, Academy of Sciences of the Czech Republic (ASCR), CzechiaPhatu William Mashela, University of Limpopo, South Africa
Copyright © 2023 Palevich, Maclean and Scott. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Nikola Palevich, nik.palevich@agresearch.co.nz
†These authors have contributed equally to this work