Transcriptional patterns of sexual dimorphism and in host developmental programs in the model parasitic nematode Heligmosomoides bakeri

Abstract

Background

Heligmosomoides bakeri (often mistaken for Heligmosomoides polygyrus) is a promising model for parasitic nematodes with the key advantage of being amenable to study and manipulation within a controlled laboratory environment. While draft genome sequences are available for this worm, which allow for comparative genomic analyses between nematodes, there is a notable lack of information on its gene expression.

Methods

We generated biologically replicated RNA-seq datasets from samples taken throughout the parasitic life of H. bakeri. RNA from tissue-dwelling and lumen-dwelling worms, collected under a dissection microscope, was sequenced on an Illumina platform.

Results

We find extensive transcriptional sexual dimorphism throughout the fourth larval and adult stages of this parasite and identify alternative splicing, glycosylation, and ubiquitination as particularly important processes for establishing and/or maintaining sex-specific gene expression in this species. We find sex-linked differences in transcription related to aging and oxidative and osmotic stress responses. We observe a starvation-like signature among transcripts whose expression is consistently upregulated in males, which may reflect a higher energy expenditure by male worms. We detect evidence of increased importance for anaerobic respiration among the adult worms, which coincides with the parasite’s migration into the physiologically hypoxic environment of the intestinal lumen. Furthermore, we hypothesize that oxygen concentration may be an important driver of the worms encysting in the intestinal mucosa as larvae, which not only fully exposes the worms to their host’s immune system but also shapes many of the interactions between the host and parasite. We find stage- and sex-specific variation in the expression of immunomodulatory genes and in anthelmintic targets.

Conclusions

We examine how different the male and female worms are at the molecular level and describe major developmental events that occur in the worm, which extend our understanding of the interactions between this parasite and its host. In addition to generating new hypotheses for follow-up experiments into the worm’s behavior, physiology, and metabolism, our datasets enable future more in-depth comparisons between nematodes to better define the utility of H. bakeri as a model for parasitic nematodes in general.

Graphical Abstract

Background

Parasitic nematodes (roundworms) infect over one billion people, causing significant morbidity [1]. They also pose a significant threat to livestock, costing billions of dollars annually in production losses and treatment costs [2, 3]. Yet, nematodes that parasitize humans or livestock are challenging to study in a controlled laboratory environment because the host is required to complete the parasite life cycle. Studying the in vivo behavior of the parasite during its infection requires samples to be taken from an infected host, not just obtaining eggs from feces. These samples can be technically challenging and expensive to obtain from large animals and not possible to obtain from humans. Heligmosomoides bakeri is well suited to being a laboratory model because it is a natural parasite of mice, which are easily maintained in a laboratory environment, and there are protocols for maintaining H. bakeri in a laboratory environment. It is closely related to economically important parasites of livestock and the hookworm parasites of humans [4, 5]. Furthermore, it can establish chronic infections in its mouse host, coarsely mimicking many other parasitic infections in a variety of hosts. Unravelling the responses to infection, adaptations to parasitism, and immunomodulatory strategies of H. bakeri will enable further development of this worm as a model for parasitic nematodes. Increased understanding of the biology of H. bakeri and what can be generalized to other related nematodes will then also enable preliminary experiments in H. bakeri to inform targeted experiments in harder to study parasites of large animals.

Much of the literature on H. bakeri refers to the species as Heligmosomoides polygyrus. While there has historically been confusion on which name to use (a third name, Nematospiroides dubius has also been used), molecular studies demonstrate that the isolates widely used in laboratories are sufficiently genetically different from the wild populations to be designated as a separate species [6,7,8]. We, therefore, refer to the laboratory strain of this parasite as H. bakeri.

Heligmosomoides bakeri has a direct life cycle with free-living and parasitic stages (Fig. 1). Eggs present in feces hatch after 36–37 h, releasing the L1 larvae, which molt to L2 larvae 28–29 h after hatching. After another 17–20 h, the larvae partially molt resulting in ensheathed L3 worms that are the infective stage of the parasite. Once eaten by a rodent host, the L3 larvae exsheath and within 24 h have penetrated the intestinal mucosa. The L3 larvae molt 90–96 h (~ 4 days) after infection into L4 larvae, which continue to develop in a granuloma formed by the host. A final molt 144–166 h (~ 7 days) after infection results in adults which migrate out to the intestinal lumen [9]. The adult worms coil around the villi where they feed, mate, and lay the eggs, which get passed with the feces [10]. The adults reside in the small intestine until they are expelled by the mouse, which can take anywhere from four to greater than 20 weeks depending on the strain of mouse [11].

Fig. 1

Life cycle of H. bakeri. Eggs present in mouse feces develop through two larval stages and arrest during their infective L3 stage. Upon being eaten by a mouse host the L3s exsheath in the stomach and upper duodenum and burrow into the duodenum tissue to form a granuloma. Development continues until the final molt to the adult form when the worms migrate to the lumen of the duodenum where they mate and lay their eggs, which get passed with the feces

Most studies on H. bakeri are motivated to understand its effects on its host’s immune system; comparatively little is known about the worm itself. Two draft genome assemblies [12, 13] are available for H. bakeri, which not only enable comparative genomic studies, but also allow for explorations of gene expression. A small set of 52 genes in H. bakeri (including chitinase, lysozyme, and glutathione S-transferases) was found to be differentially expressed between germ- and pathogen-free mice using bulk RNA-seq of mixed-sex cultures, suggesting sensing of the microbial environment by the adult worms [14]. However, no transcriptomic information exists for other stages throughout the life cycle or for sex-separated adults. Finally, attempts have been made to identify excreted/secreted products from the worm that may have immunomodulatory potential, such as miRNAs within exosomes [15] and proteins from separated and mixed-sex cultures [16].

Here, we have generated biologically replicated mRNA transcriptomes for both male and female H. bakeri at four time points throughout their infection using RNA-sequencing (RNA-seq). The samples span the fourth larval and adult stages of the worm (5, 7, 10, and 21 days post-infection). Our dataset allows us to examine, for the first time, sexual dimorphism in the L4s and adults of this species at the level of gene transcription. It also enables us to describe the major developmental changes and processes that occur in the worms during their parasitic phase. Moreover, this dataset provides the first publicly available resource to query the overall level of transcription of any gene of interest in this species throughout the parasitic stages.

Materials and methods

Mice and parasites

Male C57Bl/6 mice aged 6–8 weeks (bred and maintained at the animal care facility, Department of Biological Sciences, University of Calgary) were used. All animal experiments were approved by the University of Calgary’s Life and Environmental Sciences Animal Care Committee (protocol AC17-0083). All protocols for animal use and euthanasia were in accordance with the Canadian Council for Animal Care (Canada). Infected mice were orally gavaged with 300 third-stage H. bakeri larvae (maintained in house, original stock was a gift from Dr. Lisa Reynolds, University of Victoria, Canada) and euthanized at either 5, 7, 10, or 21 days post-initial infection. Worms were removed from the intestinal tract, placed in Dulbecco’s Modified Eagle’s Medium high glucose (Sigma cat. D5796) where they were sexed and counted. The number of worms in each sample is shown in Table 1. Worms were snap frozen and kept at −80 °C until RNA isolation.

Table 1 Metadata of the RNA-seq samples in this study

RNA isolation and quality control

Worms were lysed by adding 100 µl Trizol and homogenizing in dry ice three times, for a total of three freeze-thaw cycles and a final volume of 300 µl. The mixture was centrifuged, and the supernatant was collected into a new tube where RNA was extracted using the Zymo Direct-zol RNA Miniprep kit (cat. no. R2050) following the manufacturer’s instructions. RNA was digested once with DNase during the isolation, which has been found to sufficiently deplete DNA so as to not interfere with downstream RNA-seq analyses [17]. The RNA was then cleaned using the Zymo RNA Clean and Concentrator-5 kit (cat. no. R1015) following the manufacturer’s instructions. Quantity and quality of the total RNA were assessed on an Agilent 2200 TapeStation RNA ScreenTape following the manufacturer’s instructions.

Library preparation and sequencing

Libraries were made from suitable RNA samples using the NEBNext Ultra II Directional RNA Library Prep Kit following the manufacturer’s instructions. Libraries were multiplexed and paired-end sequenced on an Illumina NovaSeq with a S2 flow cell or SP flow cell for 300 cycles (2 × 150 bp) using the v1.5 reagent kit, following the manufacturer’s instructions. The resulting reads were deposited in the SRA under accession number PRJNA750155.

Aligning RNA-seq reads and counting reads per transcript

RNA-seq reads were aligned to the H. bakeri genome assembly obtained from WormBase ParaSite (PRJEB15396) using the splice-aware aligner STAR v2.7.3a [18]. The resulting BAM alignment file was used as input for the program featureCounts v2.0.3 [19] to count the fragments overlapping each transcript. Second stranded counts were used only on read pairs with both ends aligned in the proper orientation, and fragments overlapping multiple transcripts were fractionally counted among all matches. The commands used in all analyses can be found in Additional file 4.

Differential gene expression analysis

Read counts per transcript were rounded to the nearest integer, and differential transcript expression was analyzed using DESeq2 v1.30.1 [20]. Each sample type (e.g. D10 female or D5 male) was used as a condition for the DESeq2 object model to model the effects of the multiple conditions of worm sex and age, according to the recommendations in the DESeq2 vignette (http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#interactions). Within each subsequent pairwise comparison considered, transcripts with a false discovery rate adjusted p-value of < 0.05 were considered differentially expressed.

For male/female comparisons (to yield the sets of transcripts MS—male-specific, FS—female-specific, URM—upregulated in males, and URF—upregulated in females), pairwise comparisons considered included D5M versus D5F, D7M versus D7F, D10M versus D10F, and D21M versus D21F. Transcripts that were differentially expressed in all four pairwise comparisons were selected. Those that had ≥ tenfold higher expression in the male samples of the pairwise comparisons were retained as the MS transcripts. Those that had ≥ tenfold higher expression in the female samples of the pairwise comparisons were retained as the FS transcripts. Those that had higher expression in the male samples of the pairwise comparisons and were in unsigned network modules (see next section on co-expression network analysis) M3 or M5 and were in signed network modules M3 or M4 were retained as the URM transcripts. Those that had higher expression in the female samples of the pairwise comparisons and were in unsigned network modules M1 or M9 and were in signed network modules M1, M6, or M9 were retained as the URF transcripts.

For developmental comparisons (to yield the sets of transcripts URTDM—upregulated in tissue-dwelling males, URTDF—upregulated in tissue-dwelling females, URLDM—upregulated in lumen-dwelling males, and URLDF—upregulated in lumen-dwelling females) counts per transcript in the DESeq2 object were VST-transformed to yield expression estimates for every transcript in every sample. For males and females separately, the expression estimates were then compared between every replicate of D5 and D7 samples versus every replicate of D10 and D21 samples to select all transcripts that were upregulated in the tissue-dwelling phase or the lumen-dwelling phase. Transcripts were then filtered using the co-expression modules (URTDM—unsigned modules M2, M3, or M6 and signed modules M2, M6, M7, or M8; URLDM—unsigned modules M3, M4, or M5 and signed modules M3, M4, or M5; URTDF—unsigned modules M2, M3, M5, M6, M8, or M10 and signed modules M2, M3, M4, M7, M8, or M10; URLDF—unsigned modules M1, M3, M4, M5, or M7 and signed modules M1, M3, M5, M6, or M9). Transcripts that were differentially expressed in at least one of the D7 versus D10 or D5 versus D21 pairwise comparisons were retained.

Co-expression network analysis

Read counts per transcript were normalized into fragments per kilobase per million mapped reads (FPKM) using the R package countToFPKM v1.0. The top 75% of transcripts by expression value were then grouped into co-expression modules with the R package CEMiTool v1.14.1 [21] using the FPKM values, using the variance-stabilizing transformation and the pearson correlation. Modules were visualized in R to allow determination of the pattern of module enrichment across the sample groups.

Identification of orthologs in C. elegans, gene ontology functional enrichment analysis, Anthelmintic target identification, TPM calculation

All orthologs between H. bakeri and C. elegans were retrieved from WormBase Parasite. When transcripts of interest were identified in H. bakeri for which there was a one-to-one ortholog in C. elegans, the function of the C. elegans gene was identified through manual literature search.

GO terms enriched in any set of transcripts of interest were identified using the R package gprofiler2 v0.2.1 [22]. gprofiler2, as implemented in R, performs enrichment analysis at the gene level, even if the query set contains genes that have multiple transcripts. Consequently, when sets of transcripts of interest contained multiple isoforms of the same locus, that locus was only counted once during the enrichment analysis.

Orthologs of known anthelmintic target genes, as previously identified in C. elegans [23], were retrieved from WormBase Parasite. Additional targets of interest were taken from [24]. With the H. bakeri locus tags, expression of each target was examined. TPMs (transcripts per million) were calculated using a bash script (See Additional file 4). For average TPMs, arithmetic means of the TPMs for all replicates within a sample group were calculated in R along with the standard error.

Bead feeding assay

Adult worms were removed from the intestinal tract 10 days post-initial infection and isolated using a modified Baerman apparatus. They were washed three times with water and placed in Dulbecco’s modified Eagle’s medium-high glucose (Sigma cat. D5796) where they were sexed and counted. Worms were used immediately. Mixed sex groups of worms were incubated with or without Fluoresbrite YG microsphere beads (108 beads per 5 ml medium, Polysciences cat. 18859-1) for two days. Worms were killed in 70% ethanol. Images were acquired using the Zeiss Axio ZoomV.16 stereoscopic microscope (Bio Core Facility, Department of Biological Sciences, University of Calgary). Bright-field and fluorescence images were taken under the PlanNeoFluar Z 1×/0.25 FWD 56-mm objective with the AxioCam High-Resolution colour (HRc) and High-Resolution mono (HRm) cameras. Beads were counted manually from the images. Any regular spherical fluorescent particle with consistent size (0.5 µM diameter) in the intestinal tract was counted as a bead. When a cluster of beads was observed, the fluorescent image was zoomed in to discrete individual beads as much as possible. For clumps that could not be resolved, the bead count was estimated based on the size of the cluster. The length of each worm was measured in ImageJ. The scale was set using the Analysis -> Set scale function with the scale bar on the image. A trace of the intestinal tract of each worm was drawn using the “Freehand” tool of ImageJ, and the length of the trace was measured using the Analysis -> Measure function. Each photo was analyzed independently by two people. After all the photographs were analyzed, both people went through them again together to confirm scores. Discrepancies were recounted and the final count obtained by consensus.

Results and discussion

Mapping of bulk RNA-seq data and differential gene expression (DGE)

Using the splice-aware aligner STAR, we mapped the RNA-seq reads to the H. bakeri genome assembly obtained from WormBase ParaSite (PRJEB15396). Among all the datasets, 93.26–95.62% of the reads uniquely mapped to the reference genome (Table 1), reflecting the high quality of the RNA datasets (see Additional file 1 for additional quality control discussion). Moreover, the high mapping rate (which does not include chimeric alignments) also indicates that this genome assembly, while highly fragmented in 23,647 scaffolds, contains the vast majority of the transcribed polyadenylated RNAs.

The sample groups were compared to each other using DESeq2 to find all transcripts that are statistically significantly differentially expressed between conditions. The numbers of transcripts up- and downregulated between age-matched females and males or between adjacent time points among females and males are shown in Fig. 2. In general, differences between males and females increase with age, whereas the biggest developmental differences are between the D7 and D10 time points for females and the D5 and D7 time points for males.

Fig. 2

Numbers of transcripts statistically significantly up- and downregulated in sample group comparisons. The genome annotation contains a total of 25,215 transcripts. Numbers are positioned and colored according to the comparison being described. Purple center: transcripts upregulated (top number beside up arrow) or downregulated (bottom number beside down arrow) in females versus males at the age to the left, green: transcripts up- or downregulated in D5 versus D7 females (left) or males (right), orange: transcripts up- or downregulated in D7 versus D10 females or males, blue: transcripts up- or downregulated in D10 versus D21 females or males, black: transcripts upregulated in all female (URF) or male (URM) samples at all ages, red: tissue-dwelling versus lumen-dwelling comparisons (up: URTDF, URTDM; down: URLDF, URLDM). The corresponding number of unique genes in each comparison is listed in Additional file 3: Table S29

Transcription level differences between male and female worms

In addition to the obvious significance to reproductive biology and nematode transmission, sexual dimorphism has been found to be one of the biggest differences among life cycle stages in a variety of nematodes [25,26,27,28]. H. bakeri does not have a sex-specific chromosome, like a Y chromosome as in many mammals or W chromosome in many birds. Rather, it has an XX/XO sex determination system like C. elegans [29, 30]; therefore, there are no male- or female-specific genes in H. bakeri, only male- or female-specific expression of certain genes. To explore this sex-specific gene expression, we looked for transcripts that were consistently upregulated in the males or females. These transcripts were statistically significantly differently transcribed (p_adj < 0.05 by DGE analysis, see “Materials and methods” section and Fig. 3) between males and females of the same age that had higher expression in all the male samples (for upregulated in males—URM) or in all the female samples (for upregulated in females—URF) and were found in relevant modules of the co-expression networks (Additional file 2: Fig. S3). Throughout all analyses, gene names used refer to C. elegans genes, with the corresponding H. bakeri locus tags from WormBase ParaSite (HPOL_XXXXXXXXXX) in Additional file 3: Table S27.

Fig. 3

Flow diagram of male/female comparisons. The intersections of pairwise comparisons between age-matched male and female samples were filtered to yield the sets of transcripts discussed in the main text

Transcripts important for males versus females

Among the 1084 URM transcripts from 984 unique genes (Additional file 3: Table S3) were orthologs of genes in C. elegans with demonstrated association with males or demonstrated roles in male development (Additional file 3: Table S4). Notable examples include the transcriptionally regulated male development and patterning genes mab-3 and mab-23; the male fate specification gene her-1; the male mating behavior associated genes eat-4, cil-7, and mapk-15; and the spermatogenesis-related genes spe-4, cpb-2, mib-1, fog-3, and cpb-1. Additionally, among the URM transcripts were orthologs of genes in C. elegans with roles in the regulation of gene expression, whose expression in this H. bakeri context point to mechanisms by which male- or female-associated gene expression patterns may be established and/or maintained. These include genes involved in alternative splicing like prmt-9 and rsp-8, glycosylation ZC250.2, and protein ubiquitination and protein folding spop-1 and cnx-1. Finally, among the URM transcripts are orthologs of genes in C. elegans with roles in starvation, including ser-6, gcy-35, pck-3, tre-2, and nemt-1. Functional enrichment analysis of the entire URM set revealed 32 enriched gene ontology (GO) terms that predominantly describe protein modification processes and phosphorylation/dephosphorylation (Additional file 3: Table S5).

Among the 478 URF transcripts from 431 unique genes (Additional file 3: Table S6) were orthologs of genes in C. elegans with demonstrated roles in hermaphrodite development or maternal processes (Additional file 3: Table S7). These include genes involved in development of the vulva soc-2 and had-1, development of the spermatheca nhr-6, the sex determination pathway fox-1, oogenesis csn-1, ovulation ipp-5, egg formation perm-4 and gna-2, egg laying bar-1, tmc-2, and sek-1, and maternal roles in embryonic development vha-7. In contrast to the males, URF transcripts contain different orthologs of C. elegans genes involved in alternative splicing ddx-15, glycosylation gale-1, ugt-64, and ZK632.4, and protein ubiquitination sli-1, ZK430.7, and cif-1, as well as orthologs of genes involved in translational regulation of gene expression eif-2Bα, K07A12.4, eif-3.E, and eif-3.D. Moreover, rather than starvation-associated genes, URF transcripts include an ortholog of a C. elegans gene involved in eating eat-20. Finally, among the URF transcripts are orthologs of genes in C. elegans involved in stress responses including osmotic stress nhr-1 and oxidative stress aak-2, mek-1, mlk-1, and trx-2. Functional enrichment analysis of the entire URF set revealed 18 enriched GO terms that feature peroxisomes, mitochondria, and redox processes (Additional file 3: Table S8).

As in C. elegans, more transcripts are upregulated in males than in females

To further scrutinize the sex-specific gene expression and compare to similar analyses in C. elegans, we defined male-specific transcripts according to the criteria of [31] as transcripts that were statistically significantly differently transcribed and that had ≥ tenfold higher expression levels in all male samples compared to female samples of the same age (and vice versa for female-specific transcripts). These criteria (significantly different expression between males and hermaphrodites and ≥ tenfold higher expression), when applied to C. elegans in a microarray study, identified 285 male-specific and 160 hermaphrodite-specific genes [31]. Here, in H. bakeri, among the 160 male-specific transcripts from 143 unique genes (Additional file 3: Table S2) are 103 transcripts with no annotation information, highlighting that many aspects of the male worms remain uncharacterized. Moreover, only two of these transcripts have an ortholog in the C. elegans genome: HPOL_0000701801, which is in a many-to-one relationship with W02D9.4 (an uncharacterized protein) and HPOL_0001987901 which is in a many-to-one relationship with F29B9.7 (an uncharacterized protein). The remaining 57 transcripts come from 50 genes that are shown in Additional file 3: Table S26. Functions of these genes, inferred from their annotation, include sperm production (HPOL_0001035601), male patterning (HPOL_0001902701), collagen synthesis and cross-linking (HPOL_0000750501, HPOL_0001117101, HPOL_0001117201, HPOL_0001117301, HPOL_0001117501), signalling cascades (HPOL_0000062301, HPOL_0000164101, HPOL_0000934501, HPOL_0000982501, HPOL_0001599701, HPOL_0001692501, HPOL_0001834101, HPOL_0001912801), and glycosylation (HPOL_0000317101). A high proportion of signalling cascade proteins has also been noted for spermatogenesis-enriched transcripts in a C. elegans microarray study [32]. Notably, only one transcript from one unique gene fit the criteria as female-specific (HPOL_0000787001), and it is predicted to encode a peptidase (Additional file 3: Table S26).

URM and URF reflect both known and unexpected differences between males and females

Finding orthologs of known male-related genes in URM and female-related (hermaphrodite in C. elegans) genes in URF suggests that the filtering criteria used were appropriate to successfully recover transcripts important for the males and females, respectively. The URM and URF sets are, therefore, likely to reflect genuine differences between the males and females and not merely artifacts. The different sets of transcripts used between the males and females for alternative splicing, glycosylation, and protein ubiquitination point to the use of these processes to establish and/or maintain the important differences in gene expression between the males and females as well as potentially generate male/female isoforms of the targets of these gene products. It is surprising, however, to see oxidative and osmotic stress signatures among the URF transcripts. The collection of the worms from their host environment into medium is a stressful process that involves a transfer to a higher oxygen environment (see Additional file 1 for additional discussion on intestinal oxygen) as a source of oxidative stress. Additionally, it is possible the medium does not perfectly match the osmotic conditions of the mouse tissue and/or lumen, thus providing a source of osmotic stress. However, it is unclear why the females would be more responsive to these stresses than the males. Male worms have been found to die faster and in greater proportion when exposed to oxidative stressors like arsenite for C. elegans [33] and peroxide for H. bakeri [34]. Perhaps more sensitive and stronger stress responses in the female worms contribute to their increased tolerance of the stressors.

It is also curious that we find a starvation-like expression signature among the URM transcripts. It has been reported in C. elegans that the pharyngeal pumping rate decreases in males during mating from 180 pumps per minute to 50 [35]. However, the D5 and D7 males are larvae that are still individually encysted in the intestinal tissue with no contact with a female they are not mating. Yet, the starvation-like expression signature is present at these stages as well, which contradicts a mating-induced reduction in food intake in the males as an explanation for this transcript signature. Additionally, starvation has been reported to inhibit mating behaviors in C. elegans males [36]. The infection conditions we used here to obtain the worms are also used to generate eggs to grow to infective L3 larvae to maintain the stock of worms in the laboratory, indicating that the adult worms do mate under these conditions. It is therefore unlikely that the male worms are truly starved of food. A bead-feeding assay showed that both male and female worms ingest fluorescent beads (Fig. 4 and Additional file 6), further indicating that the males are unlikely to be starving from a total lack of food. It is possible that the males ingest less food than the females (even after accounting for differences in body size) resulting in a net energy deficit during the parasitic phase of their life. Indeed, the highest intake of beads was found among females, with high variability (Fig. 4 and Additional file 6). Alternatively, the males may have a higher energy expenditure in general compared to the females, resulting in more pronounced liberation of energy stores and a higher drive to seek food, both of which are associated with starvation responses. In support of this, male C. elegans were found to have higher carbohydrate metabolism through the glycolytic pathway than hermaphrodites [31].

Fig. 4

Bead feeding assay of male and female D10 adults, both n = 14. A The number of beads ingested, normalized to intestinal length. Males and females were differently distributed (circle). B Selected pictures of a female and male D10 worm. All the images are available in Additional file 6

High sexual dimorphism at the transcriptional level may be a general feature of nematodes, even at larval stages

The sheer number of transcripts that are significantly differentially expressed between males and females, both in URM + URF but also D10 (68.92%) and D21 (67.90%) (Figs. 2, 5c, d), reflects a significant level of sexual dimorphism at the transcriptional level in this nematode. We additionally analyzed the RNA-seq datasets from Haemonchus contortus where the adult males and adult females were sequenced separately [28] and found 12,669 of 21,477 transcripts (58.99%) to be significantly differentially expressed between the adult males and females (using our same criteria as for H. bakeri). In C. elegans, a microarray study of adult males and hermaphrodites found 14,488/26,843 (53.97%) transcripts to be differentially expressed [31]. In all three nematodes, these proportions of transcripts that are differentially expressed between the adult sexes point to significant transcriptional sexual dimorphism in general. Our dataset additionally allows examination at the fourth larval stage, where sexual dimorphism is less prominent than in adults but is still notable (D5—33.58% and D7—55.74% transcripts differentially expressed between males and females; Fig. 2). A microarray study from C. elegans L2, L3, and L4 worms that compared hermaphrodites to masculinized hermaphrodites (tra-2(ar221ts); xol-1(y9) worms are strongly masculinized, fertile XX pseudomales) also found sexual dimorphism among the larval stages [37]. Though they identified fewer genes as sex-regulated than we have (possibly because they focused on somatic tissues while excluding germline tissues, and/or because of the lower sensitivity of microarrays compared to RNA-seq, and/or because of differences between the two species), they also found fewer genes to be sex regulated among the larvae compared to adults, as well as many more male-enriched genes than hermaphrodite-enriched genes: both trends that we see here in H. bakeri. The greater number of male-enriched genes was postulated to be a consequence of the suppressive nature of the TRA-1 master sex-regulator in C. elegans, a gene for which H. bakeri has a one-to-one ortholog (HPOL_0000251701). These commonalities between the two worms support the view that TRA-1 suppressing male developmental programs in XX worms is a common feature among nematodes, whether they are dioecious like H. bakeri or androdioecious like C. elegans. Moreover, the same study found little overlap between the differentially expressed genes detected at early versus late larval stages [37]. This suggests significant, unexplored, sexual dimorphism exists even in the morphologically indistinct larval stages of other nematodes. Given the transcriptional sexual dimorphism in these nematodes, male/female differences should be considered in all future experiments, especially in dioecious species like H. bakeri or H. contortus where males and females exist in roughly uniform proportions.

Fig. 5

Volcano plots of male/female comparisons at D5 (A), D7 (B), D10 (C), and D21 (D). The transcripts in URM (blue) or URF (green) are shown. Negative log2 fold changes denote higher expression in the male group while positive values indicate higher expression in the female group. Marker genes in red are (1) her-1 (HPOL_0000740701), (2) mab-3 (HPOL_0001902701), (3) spe-4 (HPOL_0000308001), (4) csn-1 (HPOL_0001830501), (5) fox-1 (HPOL_0001264601), and (6) perm-4 (HPOL_0000100601)

Developmental transcriptional changes

The time points throughout infection for the samples used here include when H. bakeri larvae are encysted in the intestinal tissue (D5 and D7) and when the adults have emerged into the lumen (D10 and D21). To identify the transcriptional changes that occur throughout this final phase of development, we defined, for both the males and the females, a set of transcripts upregulated in the adults (lumen dwelling) and a set of transcripts upregulated in the larvae (tissue dwelling). To be considered upregulated in tissue-dwelling males (URTDM), transcripts had to: (i) have higher expression in all D5 and D7 male samples than in all D10 and D21 male samples, (ii) be in relevant modules of the co-expression networks (Additional file 2: Fig. S3), and (iii) have statistically significantly different expression in either of the tissue-dwelling versus lumen-dwelling pairwise comparisons being considered (see “Materials and methods” section and Fig. 6). To be considered upregulated in the lumen-dwelling males (URLDM), transcripts had to: (i) have higher expression in all D10 and D21 male samples than in all D5 and D7 male samples, (ii) be in relevant modules of the co-expression networks (Additional file 2: Fig. S3), and (iii) have statistically significantly different expression in either of the tissue-dwelling versus lumen-dwelling pairwise comparisons being considered. The same filtering criteria were applied to the female samples to find the transcripts upregulated in the tissue-dwelling females (URTDF) and upregulated in the lumen-dwelling females (URLDF).

Fig. 6

Flow diagram of developmental comparisons. Transcripts with expression patterns of interest were filtered according to co-expression network module and statistically significant differences in expression in pairwise comparisons between tissue- and lumen-dwelling worms to yield the sets of transcripts discussed in the main text

Transcripts important for tissue versus lumen dwellers

Among the 1829 URTDM transcripts from 1700 unique genes (Additional file 3: Table S9) were orthologs of genes in C. elegans with demonstrated roles in muscle development lev-11; epithelial development ltd-1 and efn-2; cuticle synthesis and molting dpy-31, phy-2, tsp-15, fkb-3, and numerous cuticlins and collagens; and male tail development lon-8 and mab-7 (Additional file 3: Table S10). Additionally, among URTDM transcripts were orthologs of genes in C. elegans with roles in environmental sensing nep-2 and signal transduction rrc-1 and pde-6. Functional enrichment analysis of the entire URTDM set revealed 33 enriched GO terms that featured the plasma membrane and cuticle synthesis, including structural constituents and procollagen-proline dioxygenase activities (Additional file 3: Table S11). Within the 3264 URTDF transcripts from 3043 unique genes (Additional file 3: Table S12) were orthologs of genes in C. elegans with known roles in muscle development unc-52 and stn-1; epithelial development efn-2; nervous system development irx-1, grdn-1, mig-13, kal-1, mnr-1, and mig-1; cuticle synthesis and molting bus-19, phy-2, mlt-7, zmp-2, and numerous cuticlins and collagens; and regulation of body size and length mua-3, sma-6, and lon-8 (Additional file 3: Table S13). Additionally, among URTDF transcripts were orthologs of genes in C. elegans with roles in aerobic respiration and oxygen sensing ucr-2.3, isp-1, pdl-1; environmental sensing nep-2; and nutrient absorption sms-5. Functional enrichment analysis of the entire URTDF set revealed 67 enriched GO terms that featured oxygen binding and transport, cuticle synthesis, and nervous system development (Additional file 3: Table S14). The URTDM and URTDF sets share 1560 transcripts from 1458 unique genes (Additional file 3: Table S15). Functional enrichment of the shared upregulated in tissue-dwelling worm transcripts revealed 31 enriched GO terms that feature adhesion, cuticle and molting cycle, and oxygen binding and transport (Additional file 3: Table S16).

Among the 2986 URLDM transcripts from 2650 unique genes (Additional file 3: Table S17) was an ortholog of the C. elegans spermatogenesis-related gene, ubxn-2. Additionally, among the URLDM transcripts were orthologs of genes in C. elegans with roles in the cuticle col-36 and metabolism sucg-1 and tre-1 (Additional file 3: Table S18). Functional enrichment analysis of the entire URLDM set revealed 64 enriched GO terms that feature phosphorylation/dephosphorylation, ubiquitin ligase activity, and protein modification processes (Additional file 3: Table S19). Among the 5925 URLDF transcripts from 5132 unique genes (Additional file 3: Table S20) were orthologs of genes in C. elegans with known roles in adult gonad maintenance and development ippk-1, mys-1, tin-9.2, and evl-20; egg laying fahd-1 and rfp-1; maternal factors for embryonic development mel-32, par-4, dnc-4, and emb-4; and germline maintenance and gametogenesis clk-2, dvc-1, stau-1, and etr-1 (Additional file 3: Table S21). Additionally, among URLDF transcripts were orthologs of genes in C. elegans with roles in mRNA splicing prcc-1 and prp-19; translation eif-3.C; protein deubiquitination otub-1; heparan sulfate metabolism hst-2 and pst-2; and negative regulation of aerobic respiration blos-1 or response to reoxygenation stl-1. Functional enrichment analysis of the entire URLDF set revealed 262 enriched GO terms that feature regulation of gene expression (RNA processing, glycosylation, ubiquitination, histone modification, and translation terms) and cell cycle terms (Additional file 3: Table S22). The URLDM and URLDF sets share 1007 transcripts from 843 unique genes (Additional file 3: Table S23). Functional enrichment of the shared upregulated in lumen-dwelling worms transcripts revealed 59 enriched GO terms that feature cell cycle and ubiquitin protein catabolic processes (Additional file 3: Table S24).

URTDM, URTDF, URLDM, and URLDF reflect both known and unexpected differences between life cycle stages as well as between male versus female development

Since the tissue-dwelling worms are still developing, finding in URTDM and URTDF orthologs of genes in C. elegans with known roles in development suggests that the filtering criteria used were appropriate to successfully recover transcripts important for the tissue-dwelling worms. Likewise, findings in URLDM and URLDF orthologs of genes in C. elegans with known roles in gametogenesis and reproduction suggest that the filtering criteria used were appropriate to successfully recover transcripts important for the lumen-dwelling worms. It is, therefore, worth additional investigation into the importance of the role heparan sulfate metabolism may be playing in adult female worms and into whether the additional phosphorylation and dephosphorylation activities in the adult males are purely a consequence of spermatogenesis or whether this also reflects other critical male processes.

Oxygen concentration may be an important driver of behavior within the host among parasitic nematodes in general

It has been reported previously that nematodes are capable of anaerobic respiration (e.g. Ascaris suum [38], H. contortus [39], C. elegans [40]). Additionally, a study comparing three free-living nematodes (C. elegans, Pristionchus pacificus, and Panagrolaimus superbus) and one plant-parasitic nematode (Bursaphelenchus xylophilus) found the parasitic nematode survived anaerobic conditions much longer than the free-living ones [41], suggesting surviving low oxygen environments may be a particularly important adaptation for parasitic nematodes compared to their free-living relatives. The enriched GO terms involving oxygen binding, transport, and utilization among URTDM and URTDF transcripts (and the lack of such terms among the URLDM and URLDF transcripts) suggests that an increased role for anaerobic respiration, at least in H. bakeri, occurs during the transition from the L4 to the adult stage. This transition is when the worms migrate from the intestinal mucosa, a physiologically aerobic environment, to the intestinal lumen, a physiologically hypoxic environment [42] (see Additional file 1 for additional discussion on intestinal oxygen). The transition from the L4 to adult stage is also when the final molt occurs, along with synthesis of the final cuticle. This aligns with the abundance of collagens, cuticlins, and cuticle-related GO terms that are enriched in the URTDM and URTDF transcripts. Among the many steps involved in cuticle synthesis is the modification of certain proline residues in the collagen molecules to 4-hydroxyproline [43], a process requiring molecular oxygen [43, 44], and the process most likely reflected in the enriched GO terms involving procollagen-proline dioxygenase activity (Additional file 3: Tables S11, S14, S16). Characterized human enzymes that perform this function (prolyl 4-hydroxylases) catalyze the reaction l-proline + 2-oxoglutarate + O₂ → 4-hydroxyproline + succinate + CO₂ [44]. In H. bakeri, this occurs when the L4 worms are in a physiologically aerobic environment; all other molts and synthesized cuticles in H. bakeri occur outside the host in a fully aerobic environment. We hypothesize that the levels of oxygen in the intestinal lumen are insufficient to support this collagen modification (and/or H. bakeri is unable to scavenge the required oxygen), which suggests that cuticle synthesis is a strictly aerobic process. In support of this, prolyl 4-hydroxylases are also involved in oxygen sensing through their oxygen-dependent modification and destruction of the hypoxia inducible factor [44], a transcription factor that is induced in the host epithelial cells lining the intestinal tract [42]. If luminal oxygen is indeed insufficient to support prolyl 4-hydroxylase activity, the need for molecular oxygen to get through the final molt within the host could be a driver of the worms encysting in the intestinal mucosa, a phenomenon that fully exposes the worms to the host immune system (in contrast to the lumen which is beyond the reach of many immune effectors). Many parasitic nematodes encyst within host tissue (e.g. Cooperia punctata [45] and Trichinella spiralis [46], which both invade the intestinal mucosa during their development), migrate through aerobic host tissues (e.g. Ascaris lumbricoides [47], Necator americanus [48], and Nippostrongylus brasiliensis [49], which all migrate through circulatory and pulmonary tissues), or actively seek host blood (e.g. Haemonchus contortus [50] and Ancylostoma duodenale [48]), which would all provide a rich source of molecular oxygen for synthesis of the final cuticle regardless of where the adults ultimately reside in the host.

Sexually dimorphic expression throughout the life cycle combined with incomplete annotation information can confound comparisons between species

Heligmosomoides bakeri orthologs of the C. elegans IIS pathway [51] show marked differences in expression pattern between males and females across our time points (Fig. 7), suggesting that molecular pathways involved in aging may differ between males and females. It has been reported previously in C. elegans that transcription of certain isoforms (d/f) of the daf-16/FOXO transcription factor is affected in a sex-specific manner by TRA-1 to increase daf-16 activity [52]. Our results, however, show a notable decrease in the transcription of daf-16 (HPOL_0000379201) in H. bakeri females upon reaching adulthood, in contrast to steadily increasing transcription in the males as their parasitic life progresses (Fig. 7c). Notably, however, in the current H. bakeri genome annotation only one isoform is predicted for daf-16 (HPOL_0000379201), which corresponds to the b/c isoforms in C. elegans (data not shown), highlighting both the care that must be taken when extrapolating biological knowledge from model organisms like C. elegans to other organisms and that the annotations in H. bakeri still need work.

Fig. 7

Expression in H. bakeri of orthologs of genes in C. elegans implicated in aging. Average TPMs (transcripts per million) are plotted for females (left) and males (right) at each time point. Error bars represent standard error. The insulin receptor (daf-2), kinase (age-1), and daf-16/FOXO transcription factor (daf-16) of the IIS pathway are shown along with the master energy regulator (aak-2)

Expression of immunomodulatory genes and anthelmintic target genes

The majority of the male/female and developmental transcriptional responses described here center around the intestine, hypodermis, and gonad tissues of the worms. An estimate of tissue sizes in C. elegans places these three tissues as the largest in the worm, cumulatively accounting for more than three quarters of the worm’s tissue volume [53]. Since bulk RNA-seq pools together RNA from the whole worm, it is not surprising that the signatures of these tissues dominate our results. However, one phenotype of particular interest in H. bakeri, immunomodulation of its host, is not necessarily governed by these tissues. While some of the excreted/secreted products of H. bakeri with immunomodulatory activity were found to originate in the intestine (vesicles containing microRNAs) [15], the excreted/secreted proteins with described or implied immunomodulatory activity have no demonstrated tissue of origin in the worm. Excreted/secreted proteins in other parasitic worms have been found to originate from the uterine fluid and/or other sources, like the secretory apparatus [54]. Therefore, a possible source tissue for H. bakeri immunomodulatory proteins is excretory cells, which make up a very small proportion of whole C. elegans [53]. Comparative analyses may be able to uncover differential signal originating from rarer cell types like the excretory cells, especially if the differences between conditions are extreme. We therefore investigated the expression of loci that have been associated with immunomodulatory activity to see what comparisons or expression patterns might identify other immunomodulatory candidates. Most of the excretory/secretory proteins in H. bakeri that have been associated with immunomodulatory activity have not been identified but are members of certain classes of genes with common annotation descriptions (Additional file 3: Table S25). We therefore included all loci that matched each description. Expression patterns within these groups of genes vary considerably (Additional file 2: Fig. S5). Moreover, even among described immunomodulatory proteins there are opposing expression patterns that exclude any one pairwise comparison from being more likely to identify immunomodulatory candidates (Table 2, Fig. 8, and Additional file 1 for additional discussion on immunomodulatory gene expression). Hb-TGM-6 was identified as HPOL_0001864701 [55], which has higher expression in lumen-dwelling males than in tissue-dwelling males (as found with mixed-sex protein secretion) as well as sex-linked differences in expression (Fig. 8 and Additional file 2: Fig. S5). HbARI was identified as HPBE_0000813301 in the PRJEB1203 annotation [56], which when aligned against PRJEB15396 corresponds to HPOL_0001636401 (Fig. 9), which has higher expression in the tissue-dwelling worms than in the lumen-dwelling worms (Fig. 8 and Additional file 2: Fig. S5). These data demonstrate that the expression of these two families of immunomodulators is tightly controlled. Moreover, our findings are consistent with the parasite’s need for these immunomodulators during infection. HbARI interferes with the release of IL-33, an alarmin which initiates the immune response to nematode infection and whose release is triggered by tissue damage [56]. IL-33 was found to be produced early in murine infection with the nematode Trichuris muris, which promotes worm expulsion in resistant mice [57]. The concomitant production and release of HbARI early in H. bakeri infection, when the worms are in the tissue and causing damage, helps regulate the host response to the worm’s damage and allows the worms to persist long enough to make it into the lumen where they are not causing the same level of tissue damage. In contrast, Hb-TGM is needed later in infection, where it stimulates TGF-β signalling and interferes with regulatory immune cells [55].

Table 2 Expression of immunomodulatory genes in H. bakeri

Fig. 8

Expression of immunomodulatory genes in male and female H. bakeri. Average TPM (transcript per million) are plotted for females (left) and males (right) at each time point. Error bars represent standard error. The only identified immunomodulatory genes to date are plotted: HbARI (top row), HbBARI and HbBARI_Hom2 (middle row), and HbTGM-6 (bottom row)

Fig. 9

Schematic of the HbARI locus highlighting that the annotations in H. bakeri still need work. The two genome assemblies for H. bakeri, PRJEB1203 and PRJEB15396 (shown collectively as a black line), agree in this region with the exception of the area shown in the black box. The conflicting annotations for HbARI, HPBE_0000813301 (teal) and HPOL_0001636401 (green) are shown in their annotated positions and drawn to scale. HPOL_0001636401 is shorter at the 5′ end and consequently does not contain the signal peptide predicted to be in HPBE_0000813301. Since we know HbARI is a secreted protein because it was identified among HES, the PRJEB15396 annotation must be incorrect. As it is, HPOL_0001636401 is predicted to localize to the mitochondrion (data not shown). The coverage of all of the RNA-seq reads generated in this study in this region is shown as a line graph above the annotations. The HPBE_0000813301 annotation contains an extra 5′ exon that has no support in any of the RNA-seq datasets

It is likely that H. bakeri secretes other proteins with immunomodulatory activity. However, these are unlikely to have uniform expression throughout the worm, may very well be secreted by a handful of specialized cells, and are not specific to males, females, or a life cycle stage sampled here. Therefore, a higher resolution technique, like single-cell RNA-seq, would be needed to identify other immunomodulatory candidates based on their transcriptional expression patterns.

We also analyzed the expression profiles of known anthelmintic drug targets (Tables 3, 4). Of note, most of the drug targets have medium or low levels of transcription in both sexes at all stages sampled (TPMs range from 0 to 34,856.8, with the top 10% of transcript expression corresponding to a TPM ≥ 66.9, top 30% corresponding to a TPM ≥ 14.7, and top 50% corresponding to a TPM ≥ 2.8). It is unclear whether the medium and low expressions reflect medium/low expression throughout the whole worm or high expression in a few cells and very low expression in the major tissues of the worms. The three highly expressed transcripts (top 10%) are the targets of benzimidazoles and a neuropeptide GPCR. Expression of the majority of the drug targets varies significantly between at least two of the sampled time points in both males and females, which may indicate higher drug efficacy on worms of particular ages or stages of development.

Table 3 Expression of drug targets in male H. bakeri

Table 4 Expression of drug targets in female H. bakeri

Conclusion

By separating the males and females at each time point in our biologically replicated mRNA transcriptomes that span the parasitic phase of the H. bakeri life cycle, we have been able to examine sexual dimorphism in this species and examine development without masking important sex-specific signals. We have uncovered inconsistencies in the available genome annotations for H. bakeri, along with other limitations that demonstrate a need to continue the work of annotating available genome sequences for this organism. We have identified processes that are key in establishing and/or maintaining sex-specific gene expression in this worm, without having sex-specific genes, including alternative splicing, glycosylation, and ubiquitination. Additionally, we find stronger oxidative and osmotic stress responses among female worms, potentially accounting for their previously reported better survival of oxidative assault. We hypothesize that oxygen concentration may be an important factor in the encysting behavior of larval stage H. bakeri (to get through their last molt), a factor that could be more general among parasitic nematodes. We demonstrate the advantage of quantifying the levels of all transcripts in the worm by highlighting previously unknown transcription differences in immunomodulatory genes between worm sexes and life cycle stages. Finally, the level of transcriptional sexual dimorphism we observe in this species (as well as in H. contortus and C. elegans) highlights the need to consider male/female differences in the worms in future experiments with dioecious nematodes.