Hookworm genes encoding intestinal excreted-secreted proteins are transcriptionally upregulated in response to the host’s immune system
Erich M. Schwarz, Jason B. Noon, Jeffrey D. Chicca, Carli Garceau, Hanchen Li, Igor Antoshechkin, Vladislav Ilík, Barbora Pafčo, Amy M. Weeks, E. Jane Homan, Gary R. Ostroff, Raffi V. Aroian
TL;DR
This study identifies hookworm genes that respond to the host's immune system, potentially helping the parasite survive and offering new targets for vaccines or therapies.
Contribution
The study identifies 153 hookworm genes encoding excreted-secreted proteins that are upregulated in response to host immunity, suggesting novel virulence factors and vaccine targets.
Findings
153 excreted-secreted genes in A. ceylanicum are upregulated in normal hosts compared to immunosuppressed hosts.
Immunoregulated genes are enriched for conserved and rapidly evolving proteins like CAP, ASPR, and TIMP.
Male-biased expression is observed in 50.1% of positively immunoregulated intestinal genes.
Abstract
Hookworms are intestinal parasitic nematodes that chronically infect ~500 million people. How hookworms successfully overcome host protective mechanisms is unclear, but it may involve hookworm proteins that digest host tissues, or counteract the host’s immune system, or both. To find such proteins in the zoonotic hookworm Ancylostoma ceylanicum, we used mass spectrometry to identify 565 genes encoding excreted-secreted (ES) proteins from adults, and used RNA-seq to identify genes expressed both in young adults (12 days post-infection) and in intestinal and non-intestinal tissues dissected from mature adults (19 days post-infection), infecting hamster hosts that either had normal immune systems or were immunosuppressed by dexamethasone. In adult A. ceylanicum, we observed 1,670 and 1,196 genes with intestine- and non-intestine-biased expression, respectively. Comparing hookworm gene…
Genes, proteins, chemicals, diseases, species, mutations and cell lines named across the full text — each resolved to its canonical identifier and authoritative record.
Click any figure to enlarge with its caption.
Fig 1
Fig 2
Fig 3
Fig 4
Fig 5
Fig 6
Fig 7
Fig 8
Fig 9
Fig 10
Fig 11- —http://dx.doi.org/10.13039/100006492Division of Intramural Research, National Institute of Allergy and Infectious Diseases
- —http://dx.doi.org/10.13039/100006492Division of Intramural Research, National Institute of Allergy and Infectious Diseases
- —http://dx.doi.org/10.13039/100000071National Institute of Child Health and Human Development
- —Fulbright Foundation
- —http://dx.doi.org/10.13039/100007231Cornell University
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsParasites and Host Interactions · Invertebrate Immune Response Mechanisms · Parasite Biology and Host Interactions
Introduction
The hookworms Necator americanus, Ancylostoma duodenale, and Ancylostoma ceylanicum are parasitic nematodes that infect ~500 million human beings, sickening them and lowering their economic productivity [1–3]. Hookworms are remarkably long-lived parasites: an adult hookworm can feed off its human host for up to 18 years [4,5], and infections lasting 1–5 years are common [6]. Hookworms achieve this partly by feeding on the host with enzymes optimized to digest host proteins [7], and partly by suppressing the immune systems of their hosts, which would otherwise kill or expel them quickly [8]. Only one drug, albendazole, is commonly used and effective enough to be useful in mass drug administration against hookworm infections [9,10]. However, human hookworms may be developing genetic resistance to this drug [11–13], as has already happened with dog hookworms [14]. Moreover, even effective use of albendazole does not prevent endemic hookworms from reinfecting people [15]. An anti-hookworm vaccine would be an ideal way to suppress hookworm infections [16], but no such vaccine exists, in part because we do not know which gene products mediating host-parasite interactions should be targeted as antigens [17].
It is not yet certain how hookworms (or any other parasitic nematodes) dampen the host immune system, but many possible effectors of immunomodulation exist: secreted proteases [7,18,19]; secreted protease inhibitors [20–22]; large gene families encoding diverse secreted proteins as antigenic decoys [23,24]; mimics of mammalian immune proteins [25,26]; lipid immunomodulators [27–29]; and secreted exosomes with anti-host miRNAs [30,31]. These mechanisms are by no means exclusive; some or all of them could be working at once. Such complexity may have evolved over 350 million years, when vertebrates first colonized the land and became vulnerable to parasitic nematodes [32–34]. Ignorance of how hookworms immunosuppress their hosts makes it challenging to devise vaccines against hookworm disease. It also hinders the possible use of hookworms as sources of new biological reagents against immune-mediated diseases [35,36], which appear to become more frequent with lower rates of helminth infection [37].
Hookworms and other parasitic helminths have long been known to excrete or secrete proteins into their hosts. These excreted-secreted (ES) proteins have been hypothesized, and in some cases demonstrated, to enable immunosuppression and virulence [8,38–41]. ES proteins can be exported from cephalic/pharyngeal glands, intestine, or cuticle [39,42,43]. Genes encoding ES proteins have been identified in the major human hookworm N. americanus [44], the dog hookworm Ancylostoma caninum [45], and five strongylid parasitic nematodes closely related to hookworm: Angiostrongylus vasorum [46], Haemonchus contortus [47], Heligmosomoides bakeri [48–50], Nippostrongylus brasiliensis [51,52], and Ostertagia ostertagi [53] (Fig 1). Until recently [58,59], ES proteins had not been identified for the zoonotic human hookworm A. ceylanicum, which is an important zoonotic pathogen of humans and companion animals and which can complete its lifecycle in hamsters, making it a key hookworm species for laboratory studies [2,60–62].
Morphology and evolution of the hookworm Ancylostoma ceylanicum.(A) Female of A. ceylanicum by 12 days after infection of a hamster host; scale bar, 500 μm. At this stage of development, the hookworms are young adults that have only just begun blood feeding, with mature males and a few gravid females, and with little or no egg laying. (B) Female of A. ceylanicum by 19 days after infection of a hamster host; scale bar, 1 mm. At this stage of development, the hookworms are fully mature adults that have been blood-feeding for at least a week, have mated, and have begun extensive egg laying that can last for weeks in a hamster host. Both photographs are reproduced from Schwarz et al. [54]. (C) Evolutionary relationships of the hookworm A. ceylanicum to some other nematodes discussed in this paper. This phylogeny links A. ceylanicum to other Ancylostoma and Necator hookworms, to more distantly related strongylid parasitic nematodes, and to two well-studied free-living rhabditid nematodes (C. elegans and Pristionchus pacificus). Phylogenetic data are taken from Coghlan et al., van Megen et al., and Xie et al. [55,56,57]. A. ceylanicum is highlighted in gold; three parasitic species (A. caninum, N. americanus, and H. contortus) whose ES genes were compared to those of A. ceylanicum in this paper are highlighted in red. Hookworms form a single clade (marked with a yellow bar); they are a subset of a larger clade of strongylid nematodes (marked with a green bar) which, like hookworms, are parasitic. Strongylids are, in turn, a subset of rhabditid nematodes (marked with a blue bar); this clade encompasses both parasitic nematodes and the free-living nematodes C. elegans and P. pacificus. Notably, A. ceylanicum and other strongylid parasites are more closely related to C. elegans than to P. pacificus despite having highly divergent parasitic life cycles.
To find A. ceylanicum genes whose products may enable parasitism, we have identified A. ceylanicum ES proteins, while also using RNA-seq to identify A. ceylanicum genes whose products interact with the host either through intestine-biased expression or through upregulation in response to a functioning host immune system. For the latter set of genes, we hypothesized that the parasite transcriptionally activates genes in response to the host immune system in order to neutralize the host immune response. Previously, we found that changes of A. ceylanicum gene activity during infection in vivo are much more extensive than changes seen during simulated infection in vitro [54]. Here we build on that observation by correlating ES genes with intestinal and immunoregulated genes (Fig 2). We identify genes encoding immunoregulated intestinal ES proteins that may be important for virulence or immunosuppression, that are new targets for vaccines against hookworm, and that may aid in the development of therapies against immune-mediated diseases.
An overview of experimental and computational workflow.(A) Gene prediction. Because optimal analysis of RNA-seq and proteomic data in A. ceylanicum requires gene predictions with the highest possible accuracy, we used published A. ceylanicum genomic sequence and RNA-seq data to repredict a protein-coding gene set (v2.0). Because we observed that v2.0 omitted some genes for which we had ES proteomic data, we merged non-redundant members of v1.0 to v2.0 to produce a final set of protein-coding gene predictions (v2.1). We then annotated the v2.1 genes with Phobius-predicted signal and transmembrane sequences, Pfam protein domains, orthologies to genes in other parasitic nematodes (e.g., ES protein-coding genes) and other domains/terms so that functions overrepresented by ES protein-coding or differentially expressed genes could be identified. (B) RNA-seq analysis. From either normal or dexamethasone-treated (immunosuppressed) hamster hosts, we isolated A. ceylanicum larvae 12 days after infection and A. ceylanicum adults 19 days after infection, with three biological replicates per condition (hookworm life stage and hamster status). For 19-day adults, we dissected intestines from non-intestinal tissues (12-day larvae being too small to dissect). We did RNA-seq analysis on all replicates, mapping RNA-seq reads to our v2.1 gene set. In addition, we did RNA-seq analysis on key previously published A. ceylanicum data (adult males versus adult females) to allow male-biased and female-biased genes to be identified in our results. We heatmapped all RNA-seq data, identified one pair of replicates (adult intestine, from normal versus dexamethasone-treated hamsters) that was likely to have been label-swapped, and removed them before statistical analysis of RNA-seq gene expression differences. (C) ES protein identification. We grew two sets of A. ceylanicum in normal hamster hosts to adulthood (20 days after infection), removed them from their hosts, and collected two sets of excreted/secreted (ES) proteins. We performed mass spectrometry on the two ES collections, mapped the spectra to v2.1 genes, and identified ES protein-coding genes. (D) Analysis of integrated ES and RNA-seq results. We identified statistically significant overlaps between gene categories (ES genes, differentially expressed genes, genes encoding protein domains) and biologically interpreted these results.
Results
Improved A. ceylanicum gene predictions
Analyzing gene function depends crucially on the quality of gene predictions [63,64]. We thus repredicted A. ceylanicum protein-coding genes with BRAKER2 [65], merging these repredictions (“version 2.0” or “v2.0”) with our original gene predictions (“version 1.0” or “v1.0”) to produce a hybrid v2.1 set of 33,190 protein-coding genes. Recently, Uzoechi et al. have also repredicted A. ceylanicum genes using methods similar to ours [58]; their new gene set is also substantially better than earlier predictions (94.5% completeness), though less complete than our v2.0 or v2.1 gene sets (Table 1). Having predicted the v2.1 gene set, we annotated its predicted protein products with predicted N-terminal signal sequences [71], conserved protein domains [72,73], orthologies to genes in related nematode species such as N. americanus [74], and Gene Ontology (GO) terms describing biological and molecular functions [75,76]; annotations are listed in S1 Table.
Table 1: Protein-coding gene predictions for A. ceylanicum.
Identifying ES proteins
To collect ES proteins from adult A. ceylanicum hookworms, we isolated hookworms from hamster hosts 20 days after infection in two independent experiments, incubated the hookworms in protein-free culture medium for three days, concentrated and purified their supernatants, subjected their proteins to mass spectrometry, and mapped the resulting protein spectra to a nonredundant set of protein sequences. By this means, we identified 565 genes encoding ES proteins observed in at least one ES collection (1.7% of all genes), and 350 genes encoding ES proteins in both independent collections (1.05% of all genes; S1 and S2 Tables). Uzoechi et al. and Wong et al. have also identified genes encoding ES proteins from adult A. ceylanicum [58,59]; mapping their ES genes onto our v2.1 gene set, we find 955 ES genes from their analysis (2.9% of all genes), of which 444 are identical to ES genes from ours (1.3% of all genes; 79% of our ES genes). This overlap is 27-fold greater than chance (two-tailed Fisher test, p-value = 0; S3 Table), and shows high reproduciblity of ES genes in A. ceylanicum. Conversely, each ES gene set has unique members; they collectively have 1,076 genes (S2 Table).
Out of the 565 ES genes identified, we consider two subsets for detailed analysis here: either categories of ES genes for which we can show statistically significantly enhanced overlaps with other gene categories of interest (S3 Table), or individual instances of ES genes for which we have found specific published information in the literature which indicates that an individual ES gene may have a biologically informative function or effect.
A. ceylanicum ES genes encode possible host-parasite interaction proteins
Two-thirds of our ES genes (380) were predicted to encode classically secreted proteins with N-terminal signal sequences, 4.6-fold above genome-wide background (q = 4.5•10^-181^); the remaining one-third of ES genes (185) encode proteins that are presumably non-classically secreted (Tables 2 and S3) [79]. This mix of predominant but not universal N-terminal signals also exists in ES proteins from A. caninum, N. americanus, and H. contortus (S4-S6 Tables) [44,45,47], and may reflect nonclassical secretion of ES proteins through extracellular vesicles [80,81]. For ES genes that appear to encode non-classically secreted proteins, an alternative possibility remains that these genes have 5’ exons that encode signal sequences but that have been overlooked. We have no way to completely exclude this possibility without extensive manual curation of their gene structures [82,83]. However, there are two lines of evidence that incline us to think that many, if not all, predictions of non-classical ES protein secretion are correct. First, the gene predictions we used in this study were based on extensive RNA-seq data sets [54,67,68] which drove BUSCO completeness from 87.8% to 95.1% (Table 1), and thus are likely to have detected most 5’ exons. Second, noticeable fractions of ES protein-coding genes without known N-terminal signal sequences are observable not merely in our own data for A. ceylanicum (33% of 565 ES genes without a predicted signal sequence), but also for ES protein-coding genes in three other strongylids: A. caninum (63% of 315 ES genes), N. americanus (53% of 197 ES genes), and H. contortus (40% of 844 ES genes; all three ES data sets are analyzed in S5 Table). So, if the apparent non-classical secretion of ES proteins is an artifact due to exon mispredictions, this artifact is pervasive not only for our own protein-coding gene predictions of A. ceylanicum, but also for gene predictions in three other parasitic nematode genomes.
Table 2: Protein domains overrepresented in A. ceylanicum ES genes.
Compared to the whole A. ceylanicum proteome, the 565 ES genes disproportionately encoded several protein families with plausible roles in host-parasite interactions (Tables 2 and S3 and Fig 3) [40]. These included five types of proteases (aspartyl, astacin, cysteine, serine, and metallopeptidase) that may digest host tissues and blood [7,18,85–87], along with three types of protease inhibitors (Kunitz, TIL, and TIMP) that may protect A. ceylanicum from native or host proteases; both proteases and protease inhibitors may also counteract host immune responses [88,89]. Five ES genes encoded glutathione S-transferase, an enzyme thought to detoxify free heme generated during hemoglobin proteolysis and other toxins [40,90]. Sixteen ES genes encoded Stichodactyla helianthus toxin (ShK)-related proteins, which might suppress host T cells [91,92]; and nine ES genes encoded C-type lectin proteins, which might mimic mammalian immune proteins (enabling immunosuppression) or dissociate cells (enabling tissue invasion) [25,26,54,93,94].
Groups of ES genes encoding overrepresented protein domains.This UpSet plot [84] shows, for the set of 565 ES genes, occurrences of 16 protein domains from the Pfam database that are overrrepresented in that set (Table 2). The x-axis lists individual protein domains, with the numbers of ES genes that encode them; abbreviations for protein domains are given in Table 2. The y-axis lists categories of ES genes that encode a particular set of overrepresented protein domains, with the numbers of ES genes in each category. Dots in the matrix show each category’s encoded domain or domains. Although most ES genes only encode one overrepresented protein domain, instances encoding two domains do exist, and are shown by connecting two dots with a line. Not shown are an additional 242 ES genes that lacked any of these 16 overrepresented Pfam domains. Most ES genes encode only a single protein domain, and all domains are unique to some ES genes; for instance, 109 ES genes encode only CAP domains while 22 ES genes encode only ASPR domains. A minority of ES genes encode two domains at once: for instance, 29 ES genes encode both CAP and ASPR domains, and four ES genes encode both ShK and either Astacin or CAP domains.
One quarter of the 565 ES genes (140) encoded secreted venom-like allergen proteins with CAP or SCP/TAPS domains [24]; in hookworms, these are called Ancylostoma secreted proteins, activation-associated secreted proteins, or ASPs [24,95,96]. One ASP, neutrophil inhibitory factor in A. caninum (Acan-NIF), has been experimentally shown to block neutrophil adhesion [97,98]; another ASP, hookworm platelet inhibitor in A. caninum (Acan-HPI), has been shown to prevent blood clotting [99]. We observed ES genes encoding orthologs of both these ASPs (Acey-NIF-B and Acey-HPI) which may themselves be immunosuppressant and antithrombotic proteins (S1 Table). However, this leaves the other 138 ES ASP genes with only conjectural functions.
ES genes also encode three other multigene families of secreted proteins that are suspected to be involved in infection, although their molecular roles remain largely uncharacterized: ASP-related (ASPR), transthyretin-like (TTL), and secreted clade V proteins (SCVP). ASPRs were first observed in A. ceylanicum as a divergent subfamily of ASPs transcriptionally upregulated immediately after host infection [54], and were later observed in ES proteins from N. americanus [44]. TTL proteins were previously observed as a major component of ES proteins from H. contortus [47] and as multigene families in several parasitic nematodes [100–102]; one TTL protein of H. contortus (HcTTR) antagonizes goat IL-4 in vitro, and thus may be immunosuppressive in vivo [103]. SCVPs were observed in A. ceylanicum as a novel family of secreted proteins upregulated in young hookworm adults that is conserved in C. elegans and expanded in hookworms and other strongylids [54].
Although we have focused on analyzing our own ES gene set, we have also identified statistically enriched protein domains for the 955 adult ES genes identified by Uzoechi et al. and Wong et al. (both involving Washington University, and thus collectively labeled “WashU”) and for their 511 ES genes not found in our data set (“WashU-only”; S3 Table). Domains preferentially enriched in WashU-only ES genes include aminopeptidases, fatty acid- and retinol-binding (FAR) proteins, SXP/RAL-2 proteins, and C-type lectins (Table 3). Aminopeptidases may cooperate with the five other protease types noted above to enable digestion and immunosuppression. FAR proteins are a nematode-specific protein family that may suppress immune responses in plant or animal hosts by sequestering lipid signaling molecules [104–106]. The SXP/RAL-2 protein family includes homologs of an immunodominant hypodermal antigen found in several parasitic nematodes: Ac16, in the dog hookworm A. caninum [107]; both As16 and As14, in the roundworm Ascaris suum [108–110]; Ani s 5 in the ascaridoid Anisakis simplex [111]; and Ov-RAL-2 in the filarial parasite Onchocerca volvulus [112]. Among the C-type lectins enriched in WashU-only ES genes, we observed two genes encoding structural mimics of mammalian mannose receptor (Table 4) [54]; in mammals, this receptor binds soluble secreted molecules from the whipworm Trichuris suis and from A. suum, and these A. ceylanicum mimics might thus be immunomodulatory [113]. Conversely, domains preferentially enriched in our ES genes included CAP, ASPR, cysteine protease, astacin, Kunitz protease inhibitor, ShK, TTL, and SCVP domains (Table 3).
Table 3: Protein domains differentially overrepresented in WashU-only ES genes or our ES genes.
Table 4: Individual A. ceylanicum ES genes of interest.
Other ES genes we observed were potentially relevant to infection despite not being from overrepresented families (Table 4). One encodes a macin homolog that may protect hookworms from being infected by bacteria, and thus enable their long-term residence in the human gut [114,115]. Three ES homologs of A. caninum anticoagulant protein [116] could act with Acey-HPI to prevent clotting during blood feeding. Seven other ES genes encode possible immunosuppressants: two ES apyrases and one ES adenosine deaminase could suppress both blood clotting and inflammation by hydrolyzing extracellular ATP [117–119]; phospholipases A2 and D, like FAR proteins, could suppress lipid signals in immune responses [120]; the ES macrophage migration inhibitory factor Acey-MIF-1, previously observed to bind the human MIF receptor CD74 and promote monocyte migration, could alter cellular immune responses [121,122]; and an ES deoxyribonuclease II could hydrolyze extracellular DNA nets secreted by neutrophils to trap and kill large pathogens [123].
Because ES proteins are exposed to the human host’s immune system and are likely to promote infection, ES proteins are promising targets for an anti-hookworm vaccine. Twelve ES genes encoded proteins that have already been shown to give partial protection as vaccines in laboratory mammals, or to be homologous to such proteins (Table 4): the cysteine protease Acey-CP-1 [124]; the M13 peptidases Acey-MEP-6 and Acey-MEP-7 [125,126]; the glutathione S-transferase Acey-GST-1, homologous to A. caninum’s Acan-GST-1 [127] and N. americanus’ Na-GST-1 [128]; a homolog of the immunodominant hypodermal Ac-16 antigen of A. caninum [107]; two homologs of apyrases in Teledorsagia circumcincta and H. bakeri [129,130]; two homologs of fatty acid- and retinol-binding proteins in A. ceylanicum [131] and Onchocerca volvulus [132]; one homolog of Brugia malayi calreticulin [133]; one homolog of Onchocerca volvulus fructose-1,6-bisphosphate aldolase [134]; and one homolog to A. suum enolase [135]. In addition to these strict orthologs, three instances of astacin, ASP, and TTL proteins (all overrepresented among ES proteins) have given partial protection as vaccines: Ac-MTP-1 of A. caninum, Acan-ASP-2 of A. caninum, and HcTTR of H. contortus [136–138]. Finally, two WashU-only ES genes encoded homologs of an extracellular Argonaute (exWAGO) protein in H. bakeri that binds small RNAs, is imported by mouse cells during infection with possible immunomodulatory effects, and confers partial protection as a vaccine in mice [139]; a third WashU-only ES gene encoded the aspartic protease Acey-APR-1, homologous to the vaccine subunit Na-APR-1 from N. americanus [140]. These results show the vaccine potential of ES proteins, with hundreds of other ES proteins remaining untested as vaccine subunits.
Over half of A. ceylanicum ES genes are functionally conserved in related blood-feeding parasitic nematodes: they do not merely have orthologs, but ES gene orthologs (Tables 5, S1, and S7). Out of 565 A. ceylanicum ES genes, 269 (47.6%) had ES gene orthologs in the closely related hookworm A. caninum [45], 21-fold more than expected by chance (p = 6.82•10^-301^). For the more distantly related hookworm N. americanus [44], 141/565 ES genes (25.0%) had ES gene orthologs (18-fold over chance; p = 1.78•10^-137^); for the still more distantly related blood-feeding H. contortus [47], 180/565 ES genes (31.9%) had ES gene orthologs (8.9-fold over chance, p = 4.10•10^-120^). Considering all three species at once, 350/565 ES genes (61.9%) had some ES ortholog (12-fold over chance, p = 0). We found similarly high conservation for the 860 ES genes described by Uzoechi et al. (S7 Table), and overrepresented ES gene orthologies between A. ceylanicum and N. americanus were also observed by Uzoechi et al. [58]. These orthologies of A. ceylanicum ES genes in other species are predominantly to ES genes themselves (i.e., genes experimentally shown to encode ES proteins in these other species). When ES genes in other species are set aside, the remaining orthologies of A. ceylanicum ES genes fall to, or even below, background frequencies. For A. ceylanicum ES genes without known ES gene orthologs, orthology to A. caninum dropped to 1.6-fold over chance (p = 1.31•10^-29^); in N. americanus, to 1.1-fold below chance (p = 0.22); and in H. contortus, to 2.4-fold below chance (p = 2.05•10^-30^). Thus, the A. ceylanicum ES gene set predominantly encodes ES proteins not only in A. ceylanicum itself but also in related parasitic nematodes.
Table 5: Orthologies of ES genes of A. ceylanicum to those of related parasitic nematodes.
Identification of intestinal and immunoregulated genes
Because the hookworm intestine is an important source of ES proteins, we set out to characterize intestine- and non-intestine-biased expression of hookworm genes. We infected hamsters with A. ceylanicum hookworms, allowed infections to proceed for either 12 or 19 days, collected hookworms from hamster small intestines, extracted hookworm RNAs, and performed RNA-seq with 36.9 to 55.8 million reads per biological replicate (S8 Table). Hookworms at the 12-day infection stage were young adults that were too small for us to dissect further, so RNA-seq was done on whole 12-day hookworms. By 19 days, hookworms had grown into mature adults large enough to dissect, so we separated them into intestinal and non-intestinal tissues before RNA-seq. In addition, we hypothesized that hookworms would respond to the mammalian host’s immune system by increasing their expression of proteins that counteract the mammalian immune system. To detect hookworm genes regulated by the state of the host’s immune system, 12-day and 19-day infections were performed both in normal hamster hosts (nonDEX) and in hamsters whose immune systems had been suppressed by dexamethasone (DEX). We analyzed levels and changes of gene expression both for our RNA-seq data and for published A. ceylanicum RNA-seq data from adult male intestine, male adults, and female adults (Tables 6 and S1 and Figs 4-5). We defined a gene as significantly changing its activity between two biological conditions (e.g., intestinal versus non-intestinal 19-day tissues, or intestinal nonDEX versus intestinal DEX) if the gene changed its expression by at least 2-fold (log2FC ≤ -1 or ≥ 1) with a false discovery rate (p-value corrected for multiple hypothesis testing) of no more than 0.01 (FDR ≤ 0.01). For immunoregulated genes, it is important to note that upregulation of genes in a normal host is equivalent to downregulation of genes in a dexamethasone-immunosuppressed host; we have defined such genes as positively immunoregulated (and their reverse as negatively immunoregulated), but their regulation was observed by virtue of downregulation with dexamethasone.
Table 6: Numbers of A. ceylanicum genes showing differential RNA expression (significant changes of gene expression between biological conditions).
Gene expression in A. ceylanicum.Gene activity is shown for 8,130 A. ceylanicum genes with significant differences in expression between biological conditions, with biological replicates used for differential gene expression analysis on the x-axis and individual genes on the y-axis. Expression levels are in TPM (log10). Genes sharing similar patterns of expression are split into 10 clusters. Biological replicates are: young 12-day hookworm adults from normal hosts (YA_12dpi_noDEX); young 12-day adults from immunosuppressed hosts (YA_12dpi_DEX); dissected 19-day intestines from mature hookworm adults in normal hosts (intestines_19dpi_noDEX); dissected 19-day intestines in immunosuppressed hosts (intestines_19dpi_DEX); dissected 19-day non-intestinal tissues in normal hosts (non.intestines_19dpi_noDEX); dissected 19-day non-intestinal tissues in immunosuppressed hosts (non.intestines_19dpi_DEX); adult male intestine (WashU_male_intestine); adult males (WashU_adult_male); and adult females (WashU_adult_fem). On the x-axis, biological replicates are labeled as ‘condition_X’ where X is ‘01’ through ‘04’. The adult male intestine, adult male, and adult female data in Fig 4 is a reanalysis of A. ceylanicum RNA-seq data previously published by Wei et al. [67] and Bernot et al. [68]. Additional heatmaps of gene expression for all 33,190 genes and all replicates are given in S1 Fig.
Overlaps between groups of differentially expressed A. ceylanicum genes.This UpSet plot shows the 12 most abundant overlaps of differential gene expression types (Table 6). The x-axis lists types of differential gene expression, and the numbers of genes showing them. Types are abbreviated as follows: “Female”, female-biased, determined from previously published data; “Intest”, intestine-biased, 19-day; “Male”, male-biased, determined from previously published data; “N_12d_imm”, 12-day negatively immunoregulated; “N_int_imm”, 19-day negatively immunoregulated intestinal; “N_nin_imm”, 19-day negatively immunoregulated non-intestinal; “Non_int”, non-intestine-biased, 19-day; “P_int_imm”, 19-day positively immunoregulated intestinal; and “P_nin_imm”, 19-day positively immunoregulated non-intestinal. The y-axis lists categories of genes that have a particular combination of differential gene activities, with their gene numbers; dots in the matrix show each category’s type or types of expression patterns. Although many genes have only one type of differential expression, instances with two different types do exist and are shown by connecting two dots with a line. Not shown are an additional 25,302 genes that did not show any types of differential gene expression. The most common instance of genes with two differential expression patterns is a set of 957 genes with both male-biased and positively immunoregulated intestinal expression. The next instance is a set of 225 genes with both positively immunoregulated intestinal and non-intestine-biased expression; although they are immunologically upregulated in the intestine, they are more strongly expressed in non-intestinal tissues.
By these criteria, in 19-day adults we observed 1,670 genes with intestine-biased expression and 1,196 with non-intestine-biased expression, the remaining 30,324 genes being intermediate (Table 6 and Figs 4 and 6). Of the intestine-biased genes, 723 (43.3%) were orthologous to H. contortus genes with intestine-biased expression (4.4-fold above background; p = 8.18•10^-306^; S9 Table). For comparison, when we reanalyzed published RNA-seq data for male intestine versus whole male adults [67,68] as part of the same computation, we identified 635 genes with intestine-biased expression, of which only 180 were also found in our intestine-biased gene set (S9 Table). By chance alone, one would expect to observe 32 genes shared by both intestinal sets; the observed overlap was 5.6-fold above background (p = 9.93•10^-85^) but limited (only 28% of the 635 genes were shared). The observation of over 700 genes with intestine-biased expression conserved between strongylids supports our larger gene set. The disagreement of these intestine-biased gene sets could have one or more causes: different protocols for dissection and RNA extraction; different background tissues for comparison to intestinal RNA-seq (non-intestinal here, versus whole adults in the previous data); different sexes (mixed-sex here, versus males previously); and different strains of A. ceylanicum (HY135 here, versus Indian previously), which could cause different rates of RNA-seq read mapping onto our HY135 reference genome (S10 Table) [142–144].
Differential expression of intestine-biased versus non-intestine-biased genes.This volcano plot [141] compares the expression of A. ceylanicum genes (shown as dots) in 19-day (adult) intestine versus 19-day (adult) non-intestine, with gene expression being assayed by RNA-seq of dissected intestines versus dissected non-intestinal tissues. For each gene, the x-axis shows the fold change (i.e., ratio) of its expression levels in intestine versus non-intestine; intestine-biased genes are positive and non-intestine-biased genes are negative. The y-axis shows the statistical significance (corrected for multiple hypothesis testing as a false discovery ratio, or FDR) of each gene’s fold change. Both axes have logarithmic scales (log2 for fold changes on the x-axis, -log10 for FDRs on the y-axis). In our analysis, a gene is considered to be significantly intestine-biased if its fold change is ≥ 2 with an FDR of ≤0.01, and significantly non-intestinally biased if its ratio is ≤ 0.5 with an FDR of ≤0.01; dashed lines are shown for these two criteria. Genes encoding ES proteins are plotted as red dots. Other non-ES genes that are intestine-biased are plotted as gold dots; non-intestine-biased ones are plotted as green dots; and non-biased ones are plotted as dark gray dots. ES genes are found in all expression categories, including genes with non-intestinally biased expression; this is consistent with previous observations that ES proteins can be secreted not only by intestine, but also by cephalic/pharyngeal glands or cuticle [39,42,43].
When comparing A. ceylanicum in normal (nonDEX) versus immunosuppressed (DEX) hosts, we detected substantial changes of gene expression in 19-day hookworm intestinal tissues, but almost no changes of gene expression either in 12-day young adults or in 19-day non-intestinal tissues. We observed 1,951 genes upregulated in 19-day nonDEX intestine versus 19-day DEX intestine (5.9% of all genes), along with 137 downregulated genes (0.41%; Table 6 and Fig 7). In contrast, we observed only 32 genes positively or negatively immunoregulated either in 12-day young adults or in 19-day non-intestinal tissues (Table 6 and Fig 8). One explanation for this pattern is that only 19-day intestinal tissues are directly exposed to the host’s immune system through blood feeding. Although 12-day young adult hookworms inhabit the small intestine, they have only just started feeding on blood; up to that point, they instead probably eat mucosal cells [145,146]. By 19 days, hookworm adults have been feeding on blood for up to a week; however, most of their bodily contact with this blood (and, thus, with the host’s immune system) is through the lumen of their intestine (although cephalic or pharyngeal glands may interact with their host as well).
Differential expression of positively immunoregulated intestinal versus negatively immunoregulated intestinal genes.This volcano plot compares the expression of A. ceylanicum genes assayed by RNA-seq of 19-day (adult) intestines from hookworms in normal hamster hosts (with normal immune systems) versus hookworms grown in dexamethasone-treated (immunosuppressed) hamster hosts. Aside from different conditions being compared, the x-axis, y-axis, and criteria for significance are as in Fig 6. Genes encoding ES proteins are plotted as red dots. Other non-ES genes that are positively immunoregulated in A. ceylanicum intestines are plotted as gold dots; negatively immunoregulated intestinal ones are plotted as green dots; and non-biased ones are plotted as dark gray dots. ES genes are prominently visible among genes with positively immunoregulated intestinal expression.
Differential expression of positively immunoregulated non-intestinal versus negatively immunoregulated non-intestinal genes.This volcano plot compares the expression of A. ceylanicum genes assayed by RNA-seq of 19-day (adult) non-intestinal tissues from hookworms in normal hamster hosts (with normal immune systems) versus hookworms grown in dexamethasone-treated (immunosuppressed) hamster hosts. Aside from different conditions being compared, the x-axis, y-axis, and criteria for significance are as in Fig 6. Genes encoding ES proteins are plotted as red dots. Other non-ES genes that are positively immunoregulated in A. ceylanicum non-intestinal tissues are plotted as gold dots; negatively immunoregulated non-intestinal ones are plotted as green dots; and non-biased ones are plotted as dark gray dots. Although their numbers are much smaller, a visible subset of ES genes also have positively immunoregulated non-intestinal expression.
Because nematodes lack corticosteroid receptor homologs [147], our analysis assumes that dexamethasone acts solely on the host immune system, with no direct effects on A. ceylanicum. However, when free-living larvae of the intermittently parasitic nematode Strongyloides stercoralis were treated in culture with dexamethasone by Rodpai et al., 433 S. stercoralis genes showed altered RNA expression [148]. Moreover, it has been speculated that dexamethasone administered to S. stercoralis-infected mammals may act directly on the nematodes to make them hyperinfective [149,150], perhaps through binding the S. stercoralis ecdysone receptor homolog [151]. Thus, it is possible that dexamethasone causes at least part of its transcriptional effects in A. ceylanicum by acting directly on the hookworm. That being said, some differences between S. stercoralis and A. ceylanicum may make this possibility less likely. First, the drug dosage applied to S. stercoralis by Rodpai et al. was ~ 3.6 mg of dexamethasone per 25 mg of agar, whereas the drug dosage in our experiments was 48-fold lower (3.0 mg of dexamethasone/ 1,000 mg of hamster). Second, dexamethasone has been shown through heterologous glucocorticoid-responsive transgenes to reach and affect all parts of the nematode body when fed to C. elegans [152]; in contrast, the transcriptional response of A. ceylanicum to dexamethasone treatment is almost entirely confined to its intestine. Thus, although we cannot rule out the possibility of direct transcriptional regulation by dexamethasone, we think the most likely mechanism of dexamethasone action on hookworms is by altering the immune system of their hosts.
Intestine-biased genes encode likely components of food digestion, detoxification, and absorption
Whereas ES genes overwhelmingly encoded proteins with N-terminal signal sequences alone, intestine-biased genes were modestly but significantly enriched for several predicted types of secreted or transmembrane proteins (Tables 7 and S9). These included secreted proteins, proteins predicted to have one transmembrane anchor sequence after the N-terminal signal peptide, and multi-transmembrane proteins. Intestinal genes disproportionately encoded protein families relevant to digesting, detoxifying, and absorbing food (Table 7 and Fig 9): aspartyl, cysteine, aminopeptidase, and metallopeptidase proteases [7,18,85–87]; UDP-glucoronosyl or UDP-glucosyl transferases, ecdysteroid kinase-like enzymes, and ABC transporters [90,156,157]; and major facilitator and sugar transporters [158,159]. Intestinal genes also disproportionately encoded histones; this raises the question of whether mitosis persists in the adult hookworm intestine, as has been observed or inferred for A. suum and H. bakeri [153–155]. In addition to these well-studied protein families, another notably overrepresented family was Strongylid L4 Proteins (SL4Ps), first observed in A. ceylanicum as a novel family of non-classically secreted proteins upregulated in fourth-stage (L4) larvae [54]. Analyzing previously published H. contortus intestinal RNA-seq data (S11 Table), we observed that SL4P genes were also disproportionately represented among intestine-biased H. contortus genes, along with several of the better-characterized protein families implicated in digestion, detoxification, or absorption (e.g., cysteine proteases, UDP-glucoronosyl/glucosyl transferases, and major facilitator transporters; S12 Table). Moreover, the C. elegans genes numr-1 and numr-2 encode SL4P proteins that are intestinally expressed [160]. We conclude that SL4Ps have conserved intestinal expression in hookworms and other nematodes, and that (as in C. elegans) their function might be to counteract toxic environmental stresses [161,162].
Table 7: Protein domains overrepresented in A. ceylanicum genes with intestine-biased expression.
Overlaps between groups of intestine-biased genes encoding overrepresented protein domains.This UpSet plot shows, for the set of 1,670 intestine-biased genes, occurrences of 17 Pfam domains overrrepresented in that set (Table 7). The x-axis, y-axis, and matrix are as in Fig 3, except that the genes being considered are intestine-biased rather than ES protein-coding. Abbreviations for Pfam domains are given in Table 7. Not shown are an additional 1,451 intestine-biased genes that lacked any of these 17 overrepresented Pfam domains. Although most domains are solitary, three domains that may be associated with detoxification (APH, DUF1679, and EcKL) are encoded in 13 intestinally biased genes. Another 13 intestinally biased genes encode two or three histone-related domains (Histone, H2A_C, and CBFD), which might be associated with mitotic cells in the intestine of adult parasitic nematodes [153–155].
Immunoregulated genes encode signal transduction, male-associated, host-parasite, and ES proteins
The 1,951 positively immunoregulated intestinal genes were modestly, but significantly, enriched for encoding secreted proteins, proteins with a single transmembrane anchor, or both: in other words, possible secreted or cell-surface proteins (Table 8). Overrepresented protein families included possible signal transduction components: protein kinases, protein phosphatases, and SH2 scaffold proteins [163–165]. Other overrepresented families included homologs of nematode motile sperm proteins (MSPs) [166], along with two other classes of proteins associated with MSPs: MSP fiber 2 proteins (MFP2s) [167] and DUF236 proteins [168]. Although MSPs are indeed hallmark proteins of nematode sperm cells and are usually assumed to be entirely specific to sperm, there are in fact MSPs expressed in somatic cells and required for their mobility. Two different MSP homologs in Caenorhabditis elegans have cell and axonal migration RNAi phenotypes in male linker cells and hermaphroditic neurons, with one of them being expressed in linker cells [169,170]. Moreover, msp genes are expressed in C. elegans ADL chemosensory neurons [171], and are evolutionarily retained even in parthenogenetic nematode species lacking sperm [172]. Thus, nematodes express some msp genes somatically, which fits our observation here of msp, mfp2, and duf236 gene expression in dissected intestines of adult A. ceylanicum. In addition, positively immunoregulated intestinal genes disproportionately encoded protein families with possible roles in host-parasite interactions. These families included: astacin, leishmanolysin, metallopeptidase, and trypsin-like proteases; TIMP protease inhibitors; ShK-related proteins; ASPs and ASPRs; SCVPs; and immunoglobulin domain-containing proteins (Table 8). While most of these families were also enriched among ES genes, leishmanolysins, trypsins, and immunoglobulin domains were uniquely enriched here. Despite being expressed intestinally, only one of these immunoregulated genes had intestine-biased expression; instead, 258 were non-intestine-biased (3.67-fold over chance; p = 5.51•10^-78^; S13 Table).
Table 8: Protein domains overrepresented in positively immunoregulated intestinal A. ceylanicum genes.
Genes with types of immunoregulation other than positive intestinal had fewer notable protein families (S13 Table). The 137 negatively immunoregulated intestinal genes disproportionately encoded cysteine proteases and saposins, with saposins being uniquely enriched in this gene set (36-fold over chance; q = 9.94•10^-3^). At least four C. elegans saposins have antimicrobial activity in vitro [173–176], as does a saposin from the strongylid parasite T. circumcincta [177]; conversely, one saposin in A. ceylanicum has no antibacterial activity but does lyse blood cells in vitro [178]. Negatively immunoregulated hookworm saposins might have either or both functions. Of the 26 positively immunoregulated non-intestinal genes, 15 encoded ASP proteins with CAP domains (46-fold over chance; q = 7.61•10^-19^); 10 had non-intestine-biased expression (10.7-fold over chance; p = 1.12•10^-8^). Both gene sets disproportionately encoded secreted proteins; neither set had prominently sex-biased expression.
Immunoregulated genes have significantly male-biased expression
Multicellular parasites such as hookworms are generally studied to understand mechanisms of infection common to both sexes, and indeed we collected both ES proteins and RNA from mixed-sex populations. However, males of the African tick Rhipicephalus appendiculatus excrete immunoglobulin-binding proteins that the males themselves do not seem to need, and that are instead required for coinfecting female ticks to feed on host blood efficiently [179]. This instance of sexual cooperation in a parasite made us wonder whether there existed male-biased or female-biased genes among our immunoregulated gene set. Using published RNA-seq data from A. ceylanicum adult males and females [68], we observed that 3,135 and 1,547 of our 33,190 predicted A. ceylanicum genes had male- or female-biased expression, respectively (Table 6 and Fig 5). Going on to check for overlaps of these gene sets with our immunoregulated gene set, we found that 50.1% of our positively immunoregulated intestinal genes (977/1,951; 5.3-fold over background; p = 0) were also male-biased, while only 11.6% of them (227/1,951; 2.5-fold over background; p = 3.34•10^-38^) were female-biased (Table 9 and S13 and Fig 10).
Table 9: Gene categories overrepresented in A. ceylanicum 19-day immunoregulated genes.
Overlaps between A. ceylanicum 19-day immunoregulated genes and other gene groups.This UpSet plot shows, for 2,122 genes that in some way are immunoregulated in adult (19-day) A. ceylanicum, their 16 most abundant overlaps with ES genes or other gene expression sets (Table 9). Note that because only adult immunoregulated genes are considered here, the numbers of other gene categories (e.g., 1,003 male-biased genes in this plot) are lower than for the entire A. ceylanicum genome (3,135 genes; Table 6). Types are abbreviated as follows: “ES_genes”, encoding ES proteins; “Female”, female-biased, determined from previously published data; “Intest”, intestine-biased, 19-day; “Male”, male-biased, determined from previously published data; “N_int_imm”, 19-day negatively immunoregulated intestinal; “N_nin_imm”, 19-day negatively immunoregulated non-intestinal; “Non_int”, non-intestine-biased, 19-day; “P_int_imm”, 19-day positively immunoregulated intestinal; and “P_nin_imm”, 19-day positively immunoregulated non-intestinal. The y-axis lists categories of genes that belong to a particular combination of gene types, with their gene numbers; dots in the matrix show each category’s type or types. The largest overlap has 950 genes with both male-biased and positively immunoregulated intestinal expression; the largest overlap involving ES genes has 107 genes that also show positively immunoregulated intestinal expression.
Such predominantly male-biased expression, along with the enrichment of three male-associated gene families (MSP, MFP2, and DUF236) in this gene set, raises the question of whether the changes in hookworm intestinal gene expression we observed with different host immunological backgrounds were actually due to our having harvested a higher proportion of male hookworms from normal than immunosuppressed hosts, leading to a male skew in putatively immunoregulated genes. We see two reasons why such bias is not sufficient to explain our observations. First, we only observed extensive changes of gene activity in dissected intestines from 19-day adults, while observing much smaller changes in non-intestinal tissues from those same adults (Table 6 and Figs 5, 7****, and 8). If apparent immunoregulation of genes had actually been due to overlooked male bias in the hookworms we collected from normal hosts, one would expect to see at least as much (if not more) changed gene expression in non-intestinal tissues (which contained most, if not all, of the gonadal tissue from adults) as we saw from dissected intestines. The absence of significant immunoregulation in non-intestinal tissues (or in our 12-day whole animals) is inconsistent with such an artifact. Second, positively immunoregulated intestinal genes are not a simple subset of male-biased genes, as would be expected if male selection accounted for immunological changes of gene expression: 747 of these 1,951 immunoregulated genes had no sex bias, and 227 of them had female-biased expression. We conclude that the transcriptional immunoregulation we describe here, though dominated by male-biased genes, is nevertheless real. This suggests that male and female A. ceylanicum have different responses to the host immune system, and that such differences might account for the observed male-bias from our mixed-sex A. ceylanicum 19-day intestinal RNA-seq data.
The host immune system affects ES gene expression
Protein products of ES genes are thought to affect the host immune system, but whether the host immune system affects ES genes has been unclear. Comparing the ES gene set with immunomodulated gene sets, we observed significant overlaps (Table 9 and Fig 10). Out of 1,951 positively immunoregulated intestinal genes in 19-day hookworms, 153 also encoded ES proteins (27.1% of ES genes; 4.6-fold over chance; p = 1.19•10^-59^). Of 26 positively immunoregulated non-intestinal genes, 13 were also ES genes (29-fold above chance; p = 7.47•10^-17^); such expression might reflect synthesis and secretion of ES proteins by cephalic/pharyngeal glands [42].
The 153 positively immunoregulated intestinal ES genes disproportionately encoded astacin proteases, TIMP and TIL protease inhibitors, ShK-like proteins, ASPs, ASPRs, and SCVPs (Table 10 and Fig 11). Of these genes, 69 had ES gene orthologs in A. caninum (20-fold enrichment; p = 6.00•10^-72^), 24 in N. americanus (11-fold enrichment; p = 2.35•10^-18^), and 24 in H. contortus (4.4-fold enrichment; p = 1.23•10^-9^; S14 Table). As noted above, all of the above protein families may affect immunomodulation in hookworm hosts; in addition, astacin has been shown to enable tissue invasion in vitro by A. caninum [86]. Other immunoregulated ES genes also encoded proteins with possible functions in immunoregulation, antithrombosis, or digesting host tissue (S1 Table): the ASPs Acey-NIF-B [97,98] and Acey-HPI [99], two apyrases [118,119], deoxyribonuclease II [123], hyaluronidase [180,181], and superoxide dismutase [182,183].
Table 10: Protein domains overrepresented in positively immunoregulated intestinal A. ceylanicum genes that also encode ES proteins.
Overrepresented protein domains encoded by positively immunoregulated intestinal genes that also encode ES proteins.This UpSet plot shows, for the set of 153 positively immunoregulated intestinal ES genes, occurrences of seven Pfam domains overrrepresented in that set (Table 10). The x-axis, y-axis, and matrix are as in Fig 3, except that the genes being considered are not only ES protein-coding but also positively immunoregulated in the 19-day intestine. Abbreviations for Pfam domains are given in Table 10. Not shown are an additional 53 positively immunoregulated intestinal ES genes that lacked any of these seven overrepresented Pfam domains. The domain patterns are a subset of those observed for ES genes alone (Fig 3), and represent possible virulence factors in the ES gene set.
Discussion
The ability of hookworms to suppress or evade their hosts’ immune systems and feed on their hosts’ blood and other tissues, for years, has motivated identifying hookworm genes whose products interact with the host, either by encoding ES proteins or by directly contacting host blood and tissues in the hookworm intestine [40,67]. Here we have identified ES proteins, intestine-biased genes, and immunoregulated genes in the hookworm A. ceylanicum, the last of which is unique to this study. Indeed, we hypothesized that hookworms regulate genes in response to the host immune system in order to suppress it, which motivated our immunosuppression/transcriptomic study here. Strongylid parasitic nematodes interact with the immune system despite generally not being blood-feeding [92,145,184], so blood feeding is not necessary for such interactions. However, in the case of hookworms, our findings imply that interaction with the host immune system primarily occurs through contact of the hookworm intestinal lumen with host blood.
Zoonotic A. ceylanicum productively infects both humans and other mammals, making it important to humans while also enabling laboratory studies of a true hookworm that is clinically relevant [2,62]. We have analyzed both our new gene sets and previously identified gene sets [58,67,68], to identify traits encoded either by ES genes, intestinal genes, immuoregulated genes, or all three. To enable this work, we repredicted protein-coding genes in A. ceylanicum to the highest completeness so far achieved. We have greatly expanded the number of intestine-biased genes in A. ceylanicum and found them to be extensively homologous to intestine-biased genes in H. contortus [185]. We identified a mixture of functionally suggestive and poorly understood protein families in A. ceylanicum ES proteins and intestine-biased genes, observed genes that are immunoregulated in adult intestines (where they are exposed to host blood and circulating immune factors) but not other tissues, detected an unexpected male bias in positively immunoregulated intestinal genes, and defined a positively immunoregulated subset of ES genes that have a mosaic of conservation and species-specificity.
On examination, it was possible to infer biologically coherent functions from what might seem to be a jumble of hookworm ES proteins. Immunosuppressors have been long suspected to be part of the ES protein repertoire, and we observed several protein types that could have this role. We also observed a pattern of diverse multigene families encoding ES proteins (ASPs, ASPRs, TTLs, and SCVPs); in other parasitic nematodes, such a pattern has been repeatedly seen both in ES proteins and in genetic diversity between strains or sibling species [40,102,186]. It is not clear why this pattern exists; perhaps hookworms and other parasitic nematodes use diverse ES multigene families to survive unpredictably varying host immune systems by confusing their hosts with complex and varying mixtures of protein antigens, which leads parasitic nematodes to accumulate multigene families through balancing selection [187–189]. Another possible function that consistently emerged when examining ES proteins was inhibition of blood clotting. Despite being less often discussed and less obviously virulent, antithrombotics are at least as important for parasitism as immunosuppressants. Hookworm anticoagulants were first observed in 1904 [190], and are likely to be necessary for sustained blood-feeding [119,191,192].
Although we found that ES and immunoregulated genes encode proteins already known or suspected to modulate the host immune response, we also observed two gene families overrepresented in ES and immunoregulated genes that (to our knowledge) have no previously proposed immunomodulatory functions: ASP-related (ASPR) and secreted clade V proteins (SCVP; Tables 2, 8, and 10). ASPR domains (Pfam family PF17641) are related to but distinct from the CAP domains that originally defined ASP proteins (Pfam family PF00188) [24,193]. In A. ceylanicum, ASPR genes are transcriptionally upregulated at the start of infection [54], and have been observed in N. americanus ES proteins [44]. In the most recent release (WBPS19; https://parasite.wormbase.org) of ParaSite-Wormbase [194], ASPR-encoding genes are found in eight additional strongylid species: A. caninum, A. duodenale, Angiostrongylus costaricensis, A. vasorum, Cylicocyclus nassatus, Dictyocaulus viviparus, H. bakeri, N. brasiliensis, Oesophagostomum dentatum, and Strongylus vulgaris. SCVPs (Pfam family PF17619) were first observed as a novel family of secreted proteins upregulated in young A. ceylanicum adults, with no previous related protein domains; although one to a few SCVP genes are found in free-living clade V nematodes such as in C. elegans and Pristionchus pacificus, SCVP gene families are expanded in hookworms and H. contortus [54]. The most recent release (38.0; https://www.ebi.ac.uk/interpro/entry/pfam) of Pfam [193] also shows expansions of SCVP genes in the strongylids A. caninum, C. nassatus, Haemonchus placei, O. dentatum, S. vulgaris, and Teladorsagia circumcincta. Given that they observably span multiple subdivisions of strongylids [55], ASPR and SCVP genes might be newly identified immunomodulatory factors in many or all strongylid parasites.
The strong overlap that we observed between immunoregulation and male-biased gene expression raises the question of whether hookworms have sex-specialized repertoires of genes that they express during infections [195]. There has been one striking instance (in the African tick R. appendiculatus) where males and females interact differently with their hosts and the males act to promote female success during parasitism [179]. In vitro, male and female fourth-stage larvae of H. bakeri release different ES proteins when cultured together than they do when cultured separately as male-only or female-only populations, and mixed-sex ES products from H. bakeri show more immunomodulation of dendritic cells than single-sex products [196,197]. Our current data do not indicate whether such cooperation happens during hookworm infections, let alone whether the male-biased immunoregulated genes we describe are relevant to such hypothetical cooperation. However, such a bias might ensure that both sexes are present in a successful infection, which is required for reproduction. More extensive analysis to test for sexual cooperation during infections or sex-biased gene immunoregulation should be pursued for hookworms and other parasitic nematodes.
Immunoregulated genes of A. ceylanicum have orthologs in A. caninum, N. americanus, and H. contortus. The orthologs of immunoregulated ES genes also encode ES proteins, suggesting that other parasitic nematodes may also have immunoregulated genes, and some findings indicate that they actually do. For the strongylid T. circumcincta, being embedded in host mucosa upregulates ES genes and putative immunomodulatory genes [92,184]. For the strongylid H. bakeri, infecting mouse hosts with colitis induces different ES proteins than normal hosts [198]. Finally, for the strongylid N. brasiliensis, recent transcriptomic and proteomic analyses shows genes that are positively or negatively regulated during infections of normal versus immunocompromised stat6 mutant mouse hosts, a subset of which encode ES proteins [52,199]. Thus, the immunoregulation we observe here for A. ceylanicum may be a general phenomenon found in many other parasitic nematodes.
Materials and methods
Ethics statement
All animal experiments were carried out under protocols approved by the University of Massachusetts Chan Medical School. All housing and care of laboratory animals used in this study conformed to the NIH Guide for the Care and Use of Laboratory Animals in Research (18-F22) and all requirements and all regulations issued by the USDA, including regulations implementing the Animal Welfare Act (P.L. 89–544) as amended (18-F23).
Sample procurement, preparation and storage
Cultures, infections, and collections of A. ceylanicum followed published methods [54,200]. For each A. ceylanicum infection in support of RNA-seq, we purchased Syrian golden hamsters (Mesocricetus auratus) of the HsdHan:AURA strain at 3–4 weeks of age. Hamsters were provided with food and water ad libitum. For immunosuppression experiments, we immunosuppressed half of each set of hamsters by injecting them with dexamethasone (3 mg/kg) twice per week throughout the duration of the experiment [201]; the other half were given mock injections. After the first two injections (one week), both immunocompetent and immunocompromised hamsters were infected at approximately four to five weeks of age with A. ceylanicum. Infections were allowed to progress for 12 or 19 days, after which we euthanized the hamsters and collected hookworms from small intestines dissected from the hamsters. Dissected hamster intestines were put into Hank’s Balanced Salt Solution (pre-warmed to 37°C); worms were picked quickly from the dissected tissues by hand. Collected worms were snap-frozen in liquid nitrogen and stored at −80°C until use. For A. ceylanicum infections in support of ES protein mass spectrometry, we followed a similar protocol but without injections, with two collections of A. ceylanicum at 20 days after infection, and without snap-freezing of dissected hookworms. Stages of A. ceylanicum selected for RNA-seq or proteomics are based on previously described stages of growth in golden hamsters [202]. For hookworm intestinal-specific studies, from each triplicate set of 19-day post-infection hookworms isolated from hamsters as above, we dissected both hookworm intestinal and hookworm non-intestinal tissue, and extracted RNA from the dissected tissues; for the triplicates of smaller 12-day post-infection young adult worms, such dissection was not possible, so we extracted RNA from whole 12-day young adult worms.
RNA harvesting and sequencing
Total A. ceylanicum RNA was extracted from either whole worms or from dissected tissues as in Romeo and Lemaitre [54,203]. RNA-seq was done largely as in Srinavasan et al. [204]. RNA-seq libraries were built with Illumina’s TruSeq RNA Sample Prep Kit v2 executed according to manufacturer’s instructions, using 1 µg of total RNA for each sample. RNA-seq libraries were sequenced in single-end mode with read lengths of 50 nt (for all 19-day data) or 100 nt (for all 12-day data). Newly generated A. ceylanicum RNA-seq libraries are listed in S8 Table.
Protein harvesting
After 20 days of infection, A. ceylanicum hookworms were isolated from their hamster hosts and cultured in liquid culture medium for up to 3 days at 37°C in 5% CO_2_. Liquid culture medium consisted of RPMI Medium 1640 (Gibco, Cat#11835–030) supplemented with 25 mM pH 7.2 HEPES buffer, 10 µg/ml amphotericin B (Gibco, Cat# 15290–026), and 100 U/ml penicillin/streptomycin (Gibco, Cat# 1570–063), with medium sterilized with a 0.22 µm filter. Fetal bovine serum was omitted from the medium to eliminate any added protein when quantifying ES proteins by BCA assay and to eliminate contamination for proteomics. After the first 24 hours, proteins excreted or secreted from the hookworms were collected by aspirating media only. The aspirated media were centrifuged to remove debris (1500 gs for 20 minutes, at 4°C); the resulting supernatant was concentrated 10-fold using a 3kD Amicon ultra-centrifugal filter (4000 gs for 45 minutes, at 4°C). After this first 10-fold concentration, PBS was added to match the initial volume, after which the media were reconcentrated by again being centrifuged (4000 gs for 45 minutes, at 4°C); we repeated this twice, for a total of three PBS rinses and reconcentrations. After the third PBS rinse/reconcentration, total protein was quantified using a Pierce BCA Protein Assay kit. Purified ES proteins were then frozen and stored at -80°C.
General computation
Where possible, we used mamba to install and run version-controlled software environments from bioconda [205]. For reformatting or parsing of computational results, we used Perl scripts either developed for general use or custom-coded for a given analysis. All such Perl scripts (named below with italics and the suffix “.pl”) were archived on GitHub (https://github.com/schwarzem/ems_perl). Internet sources (URLs) for other software are listed in S15 Table.
Nematode genomes, transcriptomes, coding sequences, and proteomes
For transcriptomic or proteomic analyses, published genome sequences, coding sequences, proteomes, and gene annotations of relevant nematodes were downloaded from WormBase [69] or WormBase ParaSite [194,206] (S16 Table). Published RNA-seq data of A. ceylanicum [54,67,68] and H. contortus [185] were downloaded from the European Nucleotide Archive (S11 Table). Alternative gene predictions for A. ceylanicum recently published by Uzoechi et al. [58] were obtained as a Generic Feature Format Version 3 (GFF3) annotation file (https://github.com/The-Sequence-Ontology/Specifications/blob/master/gff3.md) from Young-Jun Choi [email protected] that we have archived at https://osf.io/dxfsb. We extracted protein sequences from this GFF3 via gffread 0.12.7 [207] with the arguments ‘-g [input genome sequence FASTA file] -o/dev/null -C --sort-alpha --keep-genes -P -V -H -l 93 -y [output protein FASTA file] [input gene prediction GFF3 file]’.
Heligmosomoides species nomenclature
Recent genomic analysis has shown that the gastrointestinal parasitic nematode Heligmosomoides has two distinct species: H. bakeri and H. polygyrus [102]. This reinforced previous molecular phylogenetic observations indicating that H. bakeri is a terminal sister clade to H. polygyrus [208], and that both the morphology and the the host-specificity of H. bakeri are distinct from those of H. polygyrus [209–213]. Although laboratory strains of Heligmosomoides have often been described as H. polygyrus [214], it now appears likely that many or all of these strains have actually been H. bakeri. Thus, following previous suggestions for revised nomenclature [215], we refer exclusively to H. bakeri even when citing published work that was nominally done with H. polygyrus.
RNA-seq subsampling
We found that published A. ceylanicum RNA-seq data (S11 Table) were too extensive for gene predictions by BRAKER2 because of memory limitations. To make these data usable by BRAKER2, we selected a representative subset of them with khmer [216,217]. Because khmer requires “#/1” and “#/2” suffixes for paired-end reads, we retrofitted paired-end RNA-seq reads lacking such suffixes with retroname_fastq_reads.pl. We ran normalize-by-median.py from khmer on all data (paired and unpaired) twice, with the arguments ‘-k 31 -C 100 -M 100G’, and ‘-k 31 -C 30 -M 100G’; we then ran filter-abund.py from khmer on paired-end data with the arguments ‘--variable-coverage -C 2 [k-mer hash] [paired-end reads]’. We sorted khmer-filtered data into interleaved paired- and unpaired-end read files with paired_vs_unp_fastq.or.a.pl with the arguments ‘--r1 “#0/1” --r2 “#0/2”’.
Reprediction of protein-coding genes
To repredict protein-coding genes in A. ceylanicum, we ran braker.pl in BRAKER2 2.1.6 [65] on our published repeat-softmasked A. ceylanicum genome assembly [54], with the arguments ‘--genome [genome assembly FASTA] --prot_seq [collected proteomes FASTA] --bam [sorted mapped khmer-filtered paired-end read BAM alignment] --etpmode --softmasking --cores 48 --gff3’. BRAKER2 requires an input genome sequence to have its FASTA header lines previously stripped of comments; we did this with uncomment_FASTA_headers.pl. Running BRAKER2 also required us to install: AUGUSTUS 3.4.0 [218]; BamTools 2.5.2 [219]; cdbfasta 0.99 [220]; DIAMOND 2.0.9 [221]; GeneMark-ES/ET 4.68 [222]; and SAMtools 1.12 [223]. To guide BRAKER2 gene predictions, we used our khmer-subsampled paired-end subset of published A. ceylanicum RNA-seq data, along with predicted proteomes from the related nematodes Caenorhabditis elegans [69], H. contortus [70], and N. americanus [44]. After running BRAKER2, we renamed gene, transcript, and exon names in the GFF3 prediction file with a modified version of updateBRAKERGff.py (https://github.com/Gaius-Augustus/BRAKER/issues/416) and Perl one-line commands. To extract coding DNAs (CDS DNAs) and protein sequences from the renamed GFF3, we used gffread 0.12.7 [207] as above, with the additional argument ‘-x [output CDS DNA file]’.
This gave us a new gene set (“v2.0”) that was generally superior (in completeness as assayed by BUSCO, and nonfragmentation as assayed by count_fasta_residues.pl) to our original 2015 gene set (“v1.0”). However, we observed that some v1.0 genes encoded ES proteins (as determined by mass spectrophotometric mapping) but had no genomic overlap with v2.0 genes. To ensure no valid ES genes would be overlooked, we used BEDtools 2.30.0 [224] to identify v1.0 genes whose protein-coding exon sequences (CDSes) had no genomic overlap with v2.0 CDSes, as follows. We used get_gff3_gene_subset.pl to remove GFF3 v2.0 annotations that could not be translated by gffread, and then used Perl to extract CDS annotation lines from the GFF3 files of published v1.0 and gffread-translatable v2.0 gene predictions. We reformatted the CDS-subset files from GFF3 to BED with gff2bed in BEDOPS 2.4.41 [225]. We identified overlapping coordinates of the CDS-only BED annotation files for v1.0 and v2.0 genes with intersect from BEDtools 2.30.0 [224], using the argument ‘-loj’, and used Perl to extract a list of v1.0 genes having no CDS overlaps with v2.0 genes. Given this list of non-overlapping v1.0 genes, we used extract_fasta_subset.pl to extract their encoded CDS DNAs and proteins as subsets from the full published v1.0 CDS DNA and protein sets, and added these CDS DNA and protein subsets to our previous generated v2.0 CDS DNA and protein sequences. This gave us our final hybrid prediction set (“v2.1”) of CDS DNAs and proteins on which all further analyses in this paper were performed. The v2.1 CDS DNAs and proteome have been archived at https://osf.io/dxfsb.
We generated a v2.1 GFF3 as follows. From the published v1.0 GFF3, we extracted GFF3 annotations for these non-overlapping v1.0 genes first by running extract_parasite_GFF3_subset.pl and then by selecting ‘WormBase_imported’ annotation lines with Perl. We reformatted the resulting v1.0 non-overlapping GFF3 with AGAT 1.1.0 [226]. Likewise, we reformatted the gffread-translatable v2.0 GFF3 with AGAT. We then merged the two AGAT-reformatted GFF3s to produce a single v2.1 GFF3 with uniform formatting, which we have archived at https://osf.io/dxfsb.
Assessing quality of protein-coding genes
We determined general statistics of protein-coding gene products with count_fasta_residues.pl using the arguments ‘-e -t prot’; we obtained gene counts from maximum-isoform proteome subsets generated with get_largest_isoforms.pl using the arguments ‘-t parasite’ or ‘-t maker’. We determined and compared the completeness of predicted A. ceylanicum protein-coding gene sets with BUSCO 5.2.2 using the arguments ‘--lineage_dataset nematoda_odb10 --mode proteins’ [66], which tested a given proteome against 3,131 highly conserved single-gene orthologs in nematodes.
Gene reannotation
For the protein products of our v2.1 A. ceylanicum gene set, we predicted both N-terminal signal sequences and transmembrane alpha-helical anchors with Phobius 1.01 [71], reformatting results with tabulate_phobius_hits.pl. We predicted coiled-coil domains with Ncoils 2002.08.22 [227], reformatting results with tabulate_ncoils_x.fa.pl. We predicted low-complexity regions with PSEG 1999.06.10 [228] using the argument ‘-l’ and reformatting results with summarize_psegs.pl. We identified protein domains with Pfam 35.0 [193] database with hmmscan in HMMER 3.3.2 [229,230], using the arguments ‘--cut_ga’ to impose family-specific significance thresholds, ‘-o/dev/null’ to discard text outputs, and ‘--tblout’ to export tabular outputs; Pfam results were reformatted with pfam_hmmscan2annot.pl. We also identified protein domains with interproscan.sh in InterProScan 5.57-90.0 [73] using the arguments ‘-dp -iprlookup -goterms’, and reformatting results with tabulate_iprscan_tsv.pl. We generated Gene Ontology (GO) terms [75,76], EggNOG descriptions [231], and KEGG codes [232] with EnTAP 0.10.7-beta [233] using the argument ‘--runP’, and using selected UniProt [234] and RefSeq [235] proteome databases from highly GO-annotated model organisms (S16 Table); proteome databases were generated with makedb from DIAMOND 0.9.9 [221], itself bundled with EnTAP; annotations from EnTAP were reformatted from protein to gene annotation tables with cds2gene_EnTAP_annot.pl and cds2gene_annot.pl. We identified genes encoding possible antimicrobial peptides (AMP) by mapping predictions by Irvine et al. [115] from v1.0 to v2.1 of our gene predictions. We identified orthologies between A. ceylanicum and other nematodes with OrthoFinder 2.5.4 [74] using the arguments ‘-S diamond_ultra_sens -og’; results were reformatted with prot2gene_ofind.pl and genes2omcls.pl. For all analyses except OrthoFinder, full proteomes were used; for OrthoFinder, we used maximum-isoform proteome subsets generated with get_largest_isoforms.pl. An overall annotation table was constructed from both these and RNA-seq annotations (below) with add_tab_annots.pl.
Gene annotation for related parasitic nematodes
To make consistent comparisons of ES genes (or other categories of genes) between A. ceylanicum and other species, we annotated the protein-coding genes for three related parasites (A. caninum, H. contortus, and N. americanus) using the same methods as for A. ceylanicum above; the analyzed proteomes are listed in S16 Table.
Selection of ES gene sets from related parasitic nematodes for comparative analysis
We selected three published ES gene sets from A. caninum, H. contortus, and N. americanus for comparative analysis by the following criteria. They had to have either gene identification numbers (IDs) from public genome assemblies, or, if they did not have gene IDs from valid genomes, they needed to at least have expressed sequence tag (EST) IDs for which the original sequence data were publicly available so that they could be mapped onto genomes (and thus to modern genes with proper IDs). This requirement disqualified previously published ES gene sets for H. bakeri and N. brasiliensis. Second, the gene IDs needed to be correctly mapped from protein mass spectrometric data of one species onto genes of the same species. Remarkably, this requirement disqualified the ES genes from Angiostrongylus, which were generated from A. vasorum protein mass spectrometry data but were mapped onto genes of A. cantonensis and A. costaricensis. Third, the ES genes needed to be from parasitic nematodes that had evolved parasitism in common with A. ceylanicum hookworms. This criterion disqualified ES genes from the giant roundworm Ascaris suum or the whipworm Trichuris muris, which evolved parasitism independently of hookworms and other strongylids.
RNA-seq data
For A. ceylanicum, we generated biologically triplicated RNA-seq data for each biological condition. We also analyzed published RNA-seq data for A. ceylanicum [54,67,68] and H. contortus [185], listed in S11 Table. Before analysis, we Chastity-filtered new A. ceylanicum RNA-seq reads with quality_trim_fastq.pl using the arguments ‘-q 33 -m 50’. We then quality-filtered and adaptor-trimmed new A. ceylanicum RNA-seq reads with fastp 0.20.0 [236] using the arguments ‘--dont_overwrite --detect_adapter_for_pe --n_base_limit 0 --length_required X’, with X = 50 for new A. ceylanicum data. We quality-filtered and adaptor-trimmed previously published A. ceylanicum and H. contortus RNA-seq reads with fastp 0.20.0 using the arguments ‘--dont_overwrite --n_base_limit 0 --length_required X’ with X = 50 for A. ceylanicum.
RNA-seq expression values and significances
We quantified gene expression from RNA-seq data sets with Salmon 1.9.0, generating expression values in transcripts per million (TPM) and estimating mapped read counts per gene [237]. To prevent spurious mappings of RNA-seq reads, we used full selective alignment to a “gentrome” (a CDS DNA set, treated as a target for real mappings, combined with its genome, treated as a decoy for spurious mappings), followed by quantification using Salmon in non-alignment mode [238,239]. For Salmon’s index program, we used the arguments ‘--no-version-check --keepDuplicates -t [gentrome sequence] -d [decoy list]’; for Salmon’s quant program, we used the arguments ‘--no-version-check --seqBias --gcBias --posBias --libType A --geneMap [transcript-to-gene table]’, with ‘--unmatedReads’ used for single-end data. Results from quant.genes.sf output files were reformatted with make_salmon_TPM_slices.pl.
Heatmapping of RNA-seq data
To visualize and cluster gene expression values for RNA-seq replicates, we converted expression values of 0 TPM to empirical pseudozeros with assort_tpms.pl, with each replicate’s empirical pseudozero being defined as the lowest non-zero TPM expression value observed for all genes within that replicate. This allowed logarithmic transformation of zero expression values while defining these values in a replicate-specific way. We then converted the expression values of replicates to log_10_(TPM) scores with log10_tsv_numbers.pl. Finally, we heatmapped and dendrogram-clustered all replicates with ComplexHeatmap 2.18.0 [240] using default parameters. Heatmapping showed individual non-dexamethasone (nonDEX) and dexamethasone (DEX) intestinal RNA-seq replicates that anomalously clustered with two DEX and nonDEX intestinal RNA-seq replicates, respectively (S1 Fig); these anomalous replicates were thus not used for differential gene expression analysis.
Differential gene expression analysis
Having generated RNA-seq readcounts, we used the exactTest function of edgeR 3.36.0 [241] to compute log_2_ fold-changes (log2FC) and false discovery rate (FDR) significance values [77], to identify statistically significant changes of gene expression between pairs of biological conditions in our RNA-seq data. Dispersions were computed with edgeR’s estimateDisp function using the argument ‘robust=TRUE’. For this analysis, we used all non-anomalous replicates of our A. ceylanicum RNA-seq data (S8 Table) along with published A. ceylanicum RNA-seq data from whole fourth-stage larval (L4) males, whole adult males, whole adult females, and adult male intestine (S11 Table). To increase the statistical signal for differentially expressed genes, we used filter_minimum_readcounts.pl with the argument ‘10’ to remove genes which failed to achieve 10 mapped reads in any biological replicate before submitting readcounts to edgeR. An edgeR R script (S1 File) was generated with make_edgeR_classic_scripts.pl, and run in batch mode with R 4.1.3. In summarizing edgeR results, we defined a gene as being differentially expressed (e.g., with tissue-, sex-, or condition-biased expression) if edgeR scored it as being ≥2-fold more strongly or weakly expressed in one condition than another (either log2FC ≥ 1, or log2FC ≤ -1) with an FDR of ≤ 0.01. All other genes with intermediate expression ratios between two conditions (-1 < log2FC < 1), or with expression ratios that were determined with lower significance (FDR > 0.01) were not classified as having differential gene expression. We used edgeR’s exactTest with these criteria because it is well-adapted for analyzing small numbers of biological replicates [242].
UpSet plotting of gene set overlaps
To visualize overlaps between categories of genes, we used UpSet 2 [84,243] at the Visualization Design Lab at the University of Utah (https://vdl.sci.utah.edu/upset2), uploading our data to UpSet 2 via Multinet.app (https://upset.multinet.app).
Volcano plotting of RNA-seq data
To visualize changes of gene expression between two biological conditions in our RNA-seq data, we generated volcano plots [141] with EnhancedVolcano 1.24.0 (https://github.com/kevinblighe/EnhancedVolcano). R scripts (S2 File) were run in batch mode with R 4.4.3.
RNA-seq and differential gene expression analysis for H. contortus
These were performed with published RNA-seq data for dissected adult female intestine and whole adult female bodies (S11 Table) [185] and with predicted genes from the published chromosomal reference genome (S16 Table) [70], by the same methods as used for A. ceylanicum.
Constructing a protein set for mass spectrophotometric analysis
To combine our newly predicted A. ceylanicum proteins with alternative predictions by Coghlan et al. [55], we ran cd-hit-2d from CD-HIT 4.8.1 [244] with the arguments ‘-i [v2.1 predicted proteome] -i2 [alternative published proteome by Coghlan et al.] -o [distinct Coghlan protein sequences] -d 100 -c 1.0 -M 0 -T 1 -l 5 -s 0.0 -aL 0.0 -aS 1.0’. We then joined the v2.1 proteome with distinct Coghlan sequences to yield a nonredundant A. ceylanicum protein set for LC-MS/MS analysis.
Liquid chromatography-tandem mass spectrometry (LC-MS/MS) analysis
Peptides were analyzed with an Orbitrap Exploris 480 hybrid quadrupole-Orbitrap mass spectrometer coupled to an UltiMate 3000 RSLCnano liquid chromatography system (ThermoFisher Scientific). Peptides were loaded onto an Acclaim PepMap RSLC column (75 µm × 15 cm, 2 µm particle size, 100 Å pore size, ThermoFisher Scientific) over 15 min in 97% mobile phase A (0.1% formic acid) and 3% mobile phase B (0.1% formic acid, 80% acetonitrile) at 0.5 µL/min. Peptides were eluted at 0.3 µL/min using a linear gradient from 3% mobile phase B to 50% mobile phase B over 120 min. Peptides were electrosprayed through a nanospray emitter tip connected to the column by applying 2000 V through the ion source’s DirectJunction adapter. Full MS scans were performed at a resolution of 60,000 at 200 m/z over a range of 300–1,200 m/z with an AGC target of 300% and the maximum injection time set to ‘auto’. The top 20 most abundant precursors with a charge state of 2–6 were selected for MS/MS analysis with an isolation window of 1.4 m/z and a precursor intensity threshold of 5 × 103. A dynamic exclusion window of 20 s with a precursor mass tolerance of ± 10 ppm was used. MS/MS scans were performed using HCD fragmentation using a normalized collision energy of 30% and a resolution of 15,000 with a fixed first mass of 110 m/z. The AGC target was set to ‘standard’ with a maximum injection time of 22 ms.
Mass spectrometry data analysis
Thermo RAW files were searched against our nonredundant A. ceylanicum protein set using the SEQUEST algorithm in Proteome Discoverer 2.4 (Thermo). Data were searched for peptides with tryptic specificity with a maximum of two missed cleavages. The precursor mass tolerance was set at 10 ppm and the fragment mass tolerance was set at 0.02 Da. Search parameters included carbamidomethylation at cysteine (+57.021 Da) as a constant modification and the following dynamic modifications: oxidation at Met (+15.995 Da), acetylation at protein N-termini (+42.011 Da), Met loss at protein N-termini (-131.040 Da), Met loss+acetylation at protein N-termini (-89.030 Da). The Percolator node of Proteome Discoverer was used for PSM validation at a false discovery rate of 1%. Final results for both LC-MS/MS analyses are given in S17 Table.
Mapping ES mass spectrometry hits onto gene predictions
Our protein database for analyzing ES data included both our new v2.1 A. ceylanicum predictions and any predictions by Coghlan et al. that differed by even one amino acid residue from ours. We selected ES protein hits in S17 Table whose experimental combined q-value was ≤ 0.01; for each ES sample, we extracted their gene names as separate v2.1 and Coghlan genes, yielding four ES gene lists in all.
Because gene prediction in complex eukaryotes remains challenging [63] and because the Coghlan predictions were made on a different strain of A. ceylanicum (Indian strain, US National Parasite Collection Number 102954) than our laboratory strain used for genome sequencing and gene prediction (HY135), ES proteins whose mass spectrophotometic data fit a Coghlan prediction exclusively might do so only through slight differences in predicted gene structures or through strain-specific amino acid polymorphisms. To identify such cases where a Coghlan gene was equivalent to one of our v2.1 genes, we mapped the protein-coding exons (CDSes) of Coghlan genes onto our v2.1 gene set, as follows. We downloaded Coghlan gene annotations in GFF3 format and their corresponding A. ceylanicum genome assembly from WormBase ParaSite (S16 Table). We mapped the coordinates of Coghlan gene annotations onto our A. ceylanicum genome assembly sequence (i.e., lifted them over) with Liftoff 1.6.3 [245] using the arguments ‘-copies -polish -cds’. We extracted CDS coordinate lines from the resulting liftover GFF3 file with Perl, and did likewise for the CDS coordinates from the GFF3 of our v2.1 gene prediction set. We converted Coghlan and v2.1 CDS coordinate files from GFF3 to BED format with gff2bed in BEDOPS 2.4.41 [225]. We identified overlapping Coghlan and v2.1 CDSes with intersect in BEDtools 2.30.0 [224], using the argument ‘-loj’. We extracted columns 10 and 20 from the resulting intersection file with cut (https://www.gnu.org/software/coreutils/cut) using the argument ‘-f 10,20’, and processed the results with washu.cds_loj_umass.cds_genemap.pl using the positional arguments ‘[CDS-to-gene table] [BEDtools intersect columns 10 and 20 table]’; the CDS-to-gene map file was extracted from the FASTA headers of Coghlan and v2.1 proteomes via Perl.
We then used annot_es_v01.pl with our Coghlan-to-v2.1 gene map and our four ES gene lists to produce a unified list of ES-encoding v2.1 genes.
Mapping ES genes of Uzoechi et al. onto our gene predictions
Uzoechi et al. have recently identified ES-encoding genes in A. ceylanicum, using their own reprediction of our original 2015 gene set [58]. To compare their results to ours, we first mapped their results from ES transcripts to ES genes with map_tx_to_genes.pl. We then mapped their predicted gene set onto ours using methods similar to those for the Coghlan gene set. Because their repredictions were performed on our A. ceylanicum HY135 genome sequence rather than their A. ceylanicum Indian genome (as Coghlan predictions had been), we did not need to use Liftoff to map coordinates from one genome to the other. Subsequent mapping steps (CDS extraction with Perl; GFF to BED conversion with BEDOPS gff2bed; overlap determination with BEDtools intersect; extraction and summarizing of ES gene results with annot_es_v02.pl) were as above.
Statistical significance of overlapping gene sets
To identify significant overlaps between sets of A. ceylanicum genes, we used the Perl script motif_group_fisher.pl, which computed p-values from two-tailed Fisher tests with the Perl module Text::NSP::Measures::2D::Fisher::twotailed; for multiple hypothesis testing (e.g., comparisons to sets of protein domains) this Perl script also computed q-values via the qvalue program of the MEME 5.4.1 software suite [78].
Supporting information
S1 TableAnnotations, RNA-seq expression, and differential gene expression for the A. ceylanicum v2.1 protein-coding gene set.These annotations are provided both for all 33,190 protein coding genes (in the sheet labeled ‘Acey_v2.1 annots’), and for the following subsets of genes: encoding ES proteins (‘ES genes’, 565 genes); encoding ES proteins that were detected twice (‘ES genes detected 2x’, 350 genes); intestine-biased expression, 19-day (‘Intestine-biased’, 1,670 genes); non-intestine-biased expression, 19-day (‘Non-intestine-biased’, 1,196 genes); 19-day positively immunoregulated intestinal (‘Pos.intest.immunoreg’, 1,951 genes); 19-day positively immunoregulated intestinal and also encoding ES proteins (‘Pos.intest.immunoreg + ES’, 153 genes); 19-day negatively immunoregulated intestinal (‘Neg.intest.immunoreg’, 137 genes); 19-day positively immunoregulated non-intestinal (‘Pos.non-intest.immunoreg’, 26 genes); 19-day negatively immunoregulated non-intestinal (‘Neg.non-intest.immunoreg’, 4 genes); and 12-day negatively immunoregulated (‘Neg.12d.immunoreg’, 2 genes). In practice, users of this data will have many other different subsets of annotations that they want to extract. Thus, both here and in subsequent supplementary tables (S2-S17 Tables), the data are provided as an Excel file that the reader can filter in order to select particular subsets of data. To do this in an open Excel file, first select the top row (which will contain values that will be filtered); then select the AutoFilter option in the Data menu. For each data column, drop-down arrows should appear. These are options for filtering data in each column. For a data type/column of interest, click the drop-down arrow to open its filter menu, pick criteria for the data in that column (e.g., “Contains” a gene name or “Greater than or equal to” a given number), and click the “Apply Filter” button to get only the subset of data meeting those criteria. Multiple criteria for multiple data types can be selected (and, if they are too stringent or mutually inconsistent, may end up yielding no data). Criteria can be cleared from a given data column with the “Clear Filter” button. The data columns for S1 Table are as follows. Gene identifications and equivalencies Gene: A given predicted protein-coding gene in the A. ceylanicum genome assembly. All further data columns are pertinent to that particular gene. Gene_name: A human-readable gene name, where it exists (e.g., “Acey-CP-1” instead of simply “Acey_s0154.v2.g3234”). These names are used in the main text to discuss genes of particular interest. Mapped_v1.0_gene: Any gene or genes from our previous v1.0 predictions which overlap the v2.1 gene here by least 1 nt of coding exon sequence in the A. ceylanicum genome. Although this criterion errs on the side of sensitivity and we have made no effort to filter out short overlaps, we expect that in practice this will identify extensive exon overlaps between an earlier v1.0 gene and its reprediction in v2.1. Mapped_Coghlan_gene: Any gene or genes from the previous A. ceylanicum gene predictions by Coghlan et al. [55] which, when lifted over onto our A. ceylanicum genome assembly, overlap the v2.1 gene here by least 1 nt of coding exon sequence. Mapped_Uzoechi_gene: Any gene or genes from the recent A. ceylanicum gene repredictions by Uzoechi et al. [58] which overlap the v2.1 gene here by least 1 nt of coding exon sequence. General traits of protein products. Max_prot_size: The size of the largest predicted protein product. Prot_size: This shows the full range of sizes for all protein products from a gene’s predicted isoforms. Phobius: This denotes predictions of signal and transmembrane sequences made with Phobius 1.01 [71]. ‘SigP’ indicates a predicted signal sequence, and ‘TM’ indicates one or more transmembrane-spanning helices, with N helices indicated with ‘(Nx)’. Varying predictions from different isoforms are listed. NCoils: This shows coiled-coil domains, predicted by Ncoils 2002.08.22 [227]. Both the proportion of such sequence (ranging from 0.01 to 1.00) and the exact ratio of coiled residues to total residues are given. Proteins with no predicted coiled residues are blank. Psegs: This shows what fraction of a protein is low-complexity sequence, as detected by PSEG 1999.06.10 [228]. As with Ncoils, relative and absolute fractions of low-complexity residues are shown. Pfam: Predicted protein domains from Pfam 35.0 [193], with family-specific significance thresholds. InterPro: Predicted protein domains from InterProScan 5.57-90.0 [73]. AMP: An antimicrobial peptide (AMP) gene annotation for v1.0 of our A. ceylanicum gene predictions, predicted by Irvine et al. [115], and mapped onto our v2.1 gene predictions here. GO_Biological: Annotations from the biological subset of Gene Ontology (GO) terms [75,76], generated with EnTAP 0.10.7-beta [233]. GO_Molecular: Annotations from the molecular subset of Gene Ontology (GO) terms [75,76], generated with EnTAP 0.10.7-beta [233]. GO_Cellular: Annotations from the cellular subset of Gene Ontology (GO) terms [75,76], generated with EnTAP 0.10.7-beta [233]. EggNOG_description: EggNOG descriptions [231], generated with EnTAP 0.10.7-beta [233]. EggNOG_KEGG: KEGG codes [232], generated with EnTAP 0.10.7-beta [233]. Orthologies of protein products. OFind_Summary and OFind_Full: The results for our OrthoFinder analysis [74] of orthologies between A. ceylanicum and the related nematodes A. caninum, C. elegans, H. contortus, H. bakeri, N. americanus, N. brasiliensis, Pristionchus pacificus, and T. circumcincta. For one of these species (N. americanus) two different proteomes were included: the original predicted proteome by Tang et al. [246] labeled ‘necator_orig’, and the repredicted proteome by Logan et al. [44] labeled ‘necator_rev’. Two different views of these results are given: the summary lists taxa and gene counts, while the full results give individual gene names. ES- and gene-expression-related traits. Intest_haemonchus: Orthologous H. contortus genes (taken from OFind_Full) that have intestine-biased gene expression, as computed by our analysis of previously published H. contortus RNA-seq data (S11 Table). Non-intest_haemonchus: Orthologous H. contortus genes (taken from OFind_Full) that have non-intestine-biased gene expression, as computed by our analysis of previously published H. contortus RNA-seq data (S11 Table). ES_component: A gene whose protein product was detected in at least one of our two ES mass spectrometry experiments, E20201108-05 and E20201108-07. Note that this was used as a binary (Boolean) classification for statistical analyses of gene set overlaps; see below for other such binary classifications. ES_observations: Observations of a gene’s protein product being present in either E20201108-05, or E20201108-07, or both. Coghlan-spec_ES: In our mass spectrometry analysis, our observation of a peptide mapping specifically to a Coghlan gene-encoded protein that was distinct from our v2.1 proteins by at least one amino acid residue, but whose Coghlan gene’s coding exons could then be lifted over to the coding exons of the v2.1 gene here (Mapped_Coghlan_gene). For each such mapping, the ES observations behind it are noted (either E20201108-05, or E20201108-07, or both). Uzoechi_ES: An ES-encoding gene detected by Uzoechi et al. [58] whose coding exons mapped onto the coding exons of this v2.1 gene (in Mapped_Uzoechi_gene). Note that Uzoechi et al. separately observed both male and female ES proteins, and those specific observations are given here. Wong_ES: An ES-encoding gene detected by Wong et al. [59] whose coding exons mapped onto the coding exons of this v2.1 gene (in Mapped_Uzoechi_gene). Note that Wong et al. separately observed both L3i and adult ES proteins, and those specific observations are given here. ES_a_caninum: Orthology to an A. caninum gene encoding an A. caninum ES protein observed by Morante et al. [45], extracted from OFind_Full. ES_necator_rev: Orthology to an N. americanus gene encoding an N. americanus ES protein observed either by Logan et al. [44] or by Wong et al. [59], extracted from OFind_Full. ES_haemonchus: Orthology to an H. contortus gene encoding an H. contortus ES protein observed by Wang et al. [47], extracted from OFind_Full. Binary (Boolean) classifications of genes. Male-biased: Annotation here indicates a gene with male-biased gene expression, as defined by ≥2-fold higher expression in males versus females, with a false discovery rate (FDR) of ≤ 0.01. Female-biased: Annotation here indicates a gene with female-biased gene expression, as defined by ≥2-fold higher expression in females versus males, with a false discovery rate (FDR) of ≤ 0.01. Intestine-biased: Annotation here indicates a gene with intestine-biased gene expression, as defined by ≥2-fold higher expression in 19-day intestinal tissue from normal hosts versus 19-day non-intestinal tissue from normal hosts, with a false discovery rate (FDR) of ≤ 0.01. Non-intestine-biased: Annotation here indicates a gene with non-intestine-biased gene expression, as defined by ≥2-fold higher expression in 19-day non-intestinal tissue from normal hosts versus 19-day intestinal tissue from normal hosts, with a false discovery rate (FDR) of ≤ 0.01. Any.19d.immunoreg: Annotation here indicates a gene with either intestinal or non-intestinal and either positively or negatively immunoregulated gene expression, as defined by ≥2-fold higher or lower expression in 19-day intestinal or non-intestinal tissue from normal hosts versus 19-day corresponding tissue from dexamethasone-immunosuppressed hosts, with a false discovery rate (FDR) of ≤ 0.01. Any.19d.pos.immunoreg: Annotation here indicates a gene with either intestinal or non-intestinal positively immunoregulated gene expression, as defined by ≥2-fold higher in 19-day intestinal or non-intestinal tissue from normal hosts versus 19-day corresponding tissue from dexamethasone-immunosuppressed hosts, with a false discovery rate (FDR) of ≤ 0.01. Pos.intest.immunoreg: Annotation here indicates a gene with intestinal positively immunoregulated gene expression, as defined by ≥2-fold higher expression in 19-day intestinal tissue from normal hosts versus 19-day intestinal tissue from dexamethasone-immunosuppressed hosts, with a false discovery rate (FDR) of ≤ 0.01. Neg.intest.immunoreg: Annotation here indicates a gene with intestinal negatively immunoregulated gene expression, as defined by ≥2-fold higher expression in 19-day intestinal tissue from dexamethasone-immunosuppressed hosts versus 19-day intestinal tissue from normal hosts, with a false discovery rate (FDR) of ≤ 0.01. Pos.non-intest.immunoreg: Annotation here indicates a gene with non-intestinal positively immunoregulated gene expression, as defined by ≥2-fold higher expression in 19-day non-intestinal tissue from normal hosts versus 19-day non-intestinal tissue from dexamethasone-immunosuppressed hosts, with a false discovery rate (FDR) of ≤ 0.01. Neg.non-intest.immunoreg: Annotation here indicates a gene with non-intestinal negatively immunoregulated gene expression, as defined by ≥2-fold higher expression in 19-day non-intestinal tissue from dexamethasone-immunosuppressed hosts versus 19-day non-intestinal tissue from normal hosts, with a false discovery rate (FDR) of ≤ 0.01. Neg.12d.immunoreg: Annotation here indicates a gene with negatively immunoregulated gene expression in 12-day young adults, as defined by ≥2-fold higher expression in 12-day young adults from dexamethasone-immunosuppressed hosts versus 12-day young adults from normal hosts, with a false discovery rate (FDR) of ≤ 0.01. (Note that no genes were observed with positively immunoregulated gene expression in 12-day young adults, so no ‘Pos.12d.immunoreg’ data column was needed). WashU_any_adult_ES: Annotation here indicates a gene which was detected as encoding some sort of A. ceylanicum ES protein either by Uzoechi et al. [58] or by Wong et al. [59]. Uzoechi_male.only_ES: Annotation here indicates a gene which was detected as encoding a male-specific ES protein by Uzoechi et al. [58], annotated as solely ‘Uzoechi_male_ES’ in Uzoechi_ES. Uzoechi_female.only_ES: Annotation here indicates a gene which was detected as encoding a female-specific ES protein by Uzoechi et al. [58], annotated as solely ‘Uzoechi_female_ES’ in Uzoechi_ES. Uzoechi_both_ES: Annotation here indicates a gene which was detected as encoding an ES protein in both males and females by Uzoechi et al. [58], annotated as both ‘Uzoechi_male_ES’ and ‘Uzoechi_female_ES’ in Uzoechi_ES. Gene expression. [X].TPM: For each individual RNA-seq data set (with ‘X’ denoting the data set’s abbreviation), this gives gene expression levels in TPM, computed by Salmon 1.9.0 [237]. Keys to all abbreviations are given in S8 and S11 Tables. Biological replicates of RNA-seq samples are denoted by suffixes such as ‘_1’, ‘_2’, or ‘_3’. [X].reads: For each individual RNA-seq data set (with ‘X’ denoting the data set’s abbreviation), this gives numbers of mapped RNA-seq reads per gene, computed for individual RNA-seq data sets by Salmon 1.9.0 [237], with fractional values rounded down to integers. Differential gene expression between biological conditions. [X].vs.[Y].logFC: The fold-changes of gene expression between biological condition X and biological condition Y, expressed as log2 values, and with positive values representing greater expression in condition X. The values listed here are only those were computed to be significant using edgeR 3.36.0 [241], with multiple biological RNA-seq replicates for most conditions being compared, with all biological replicates being analyzed in a single edgeR run, and with significant results annotated for individual genes. Biological conditions of RNA-seq samples are abbreviated as shown in [X].TPM or [X].reads but without the replicate suffixes. [X].vs.[Y].FDR: The false discovery rate (FDR) for gene expression changes between biological condition X and biological condition Y, annotated for individual genes. The FDR for a given set of positive results is defined as that significance threshold which, if accepted, will lead to the entire set of positives having a collective false-positive rate no greater than the FDR; it therefore provides a way to correct for testing multiple hypotheses without rejecting excessive numbers of true positives. As with [X].vs.[Y].logFC, only changes that were computed to be significant by edgeR are listed.(XLSX)
S2 TableNumbers of genes encoding adult ES proteins of *A. ceylanicum.*Numbers are given for our observations, for recently published observations by Uzoechi et al. [58] and Wong et al. [59], and for their subsets and overlaps. Uzoechi et al. independently observed ES proteins from males and females, so a gene in their set could encode ES proteins from either or both sexes. ES genes by Uzoechi et al. or Wong et al. were mapped from their predicted gene set (Table 1) to ours before analysis. Of our observed ES genes, 14 had no corresponding gene in the Uzoechi gene set, and thus might have actually been present in the Uzoechi or Wong ES mass spectrometry data yet not mapped onto a gene.(XLSX)
S3 TableStatistical analyses of the overlaps between A. ceylanicum genes encoding protein domains (or other traits) and A. ceylanicum genes encoding ES proteins.Each category of genes (e.g., genes encoding a particular protein domain or having a particular trait) was compared for its degree of overlap to a set of ES genes, and statistically analyzed for the non-randomness of this overlap by a two-tailed Fisher test. In situtations where many categories were compared at once (for instance, an entire set of protein domains from Pfam, InterPro, or Phobius), p-values initially generated by Fisher testing were used to compute q-value significance scores which corrected for multiple hypothesis testing. Statistics are provided for the set of 565 ES genes described in this paper (labeled “UMass_ES”), for the set of 955 adult ES genes described by Uzoechi et al. [58] or Wong et al. (labeled “WashU_ES”), and for the set of 511 adult ES genes described only by Uzoechi et al. or Wong et al. but not by data in this paper (labeled “WashU-only ES”). For each ES set, spreadsheets are given for analyses of overlaps with the following traits; signal or transmembrane domain organizations predicted by Phobius; protein domains in Pfam; protein domains in InterPro; and various binary comparisons with other gene traits (“Various”). In addition, two spreadsheets are given for comparisons of the statistical significances (q-values) for overrepresented Pfam and InterPro domains, with the comparisons being between our set of 565 ES genes and the WashU-only set of 511 adult ES genes. In these comparisons, only domains with a q-value ≤ 0.1 for at least one gene set were considered, and the q-values were compared by computing their ratios for each domain. For a given domain, a q-value ratio very far from 1 indicated that that domain was much more significantly overrepresented in one gene set than in the other. Each analysis of overrepresentation in gene sets provides the following data: “Motif” denotes the specific protein domain (or other gene category) for which genes encoding it are being tested for nonrandomly high (or low) overlap with the ES gene set; “All_genes” gives the total number of v2.1 protein-coding genes within which overlaps were tested; “Motif_genes” gives the number of genes annotated with the protein domain (or other trait) being tested for overlap; “Class_genes” gives the number of ES genes being tested for overlap; “Motif.Class_overlap” gives the observed number of genes falling into both categories; “Exp_rand_overlap” gives the number of genes that would be expected to overlap purely randomly; “Enrichment” gives the ratio of observed versus random overlaps (note that this ratio can be lower than 1, and in fact can be as low as 0); “p-value” gives an initial stastistical significance for the observed overlap, computed by a two-tailed Fisher test; “q-value” gives, for cases of many comparisons at once (e.g., testing for all Pfam domains simultaneously), a significance score that conservatively corrects for multiple hypothesis testing [77,78]. Note that q-values were not computed for simple binary comparisons of gene traits listed in “Various”, which are annotated for A. ceylanicum genes in S1 Table, and defined in its table legend: Intestine-biased; Non-intestine-biased; Pos.intest.immunoreg; Neg.intest.immunoreg; Pos.non-intest.immunoreg; Neg.non-intest.immunoreg; Male-biased; Female-biased; Uzoechi_any_ES; Uzoechi_male.only_ES; Uzoechi_female.only_ES; and Uzoechi_both_ES. Also, again note that highly significant overlaps can be either higher or lower than the randomly expected genome-wide background rate.(XLSX)
S4 TableAnnotations, RNA-seq expression, and differential gene expresssion for the H. contortus protein-coding gene set.We annotated predicted protein products with N-terminal signal sequences, conserved protein domains, orthologies to genes in related nematode species, and Gene Ontology (GO) terms describing biological and molecular functions. Protein-coding genes were predicted by Doyle et al. [70]. Since the annotation methods we used for this H. contortus proteome were identical to those we used for our A. ceylanicum v2.1 proteome, almost all the data columns used here are equivalent to those in S1 Table. One data column unique to this table is “Hco_ES”, which denotes H. contortus genes previously demonstrated by Wang et al. [47] to encode ES proteins.(XLSX)
S5 TableAnnotations for protein-coding gene sets of A. caninum and *N. americanus.*We annotated predicted protein products with N-terminal signal sequences, conserved protein domains, orthologies to genes in related nematode species, and Gene Ontology (GO) terms describing biological and molecular functions. Protein-coding genes were predicted for A. caninum by Coghlan et al. [55] and for N. americanus by Logan et al. [44]. Since the annotation methods we used for these proteomes were identical to those we used for our A. ceylanicum v2.1 proteome, almost all the data columns used here are equivalent to those in S1 Table. Two data columns unique to these tables are Acan_ES and Necator_ES, which respectively denote A. caninum or N. americanus genes previously demonstrated by Morante et al. [45] or Logan et al. [44] to encode ES proteins.(XLSX)
S6 TableStatistical analyses of the overlaps between A. caninum, N. americanus, or H. contortus genes encoding protein domains (or other traits) and A. caninum, N. americanus, or H. contortus genes encoding ES proteins.ES gene annotations are taken from S4 and S5 Tables. Statistical significances of overlaps between domain/trait genes and ES genes were computed and described as in S3 Table.(XLSX)
S7 TableStatistical analyses of overlaps between pairs of A. ceylanicum gene sets, with A. ceylanicum genes encoding ES proteins versus A. ceylanicum genes having various orthologies to either ES genes or any genes in A. caninum, H. contortus, or N. americanus.“ES_a_caninum”, “ES_haemonchus”, and “ES_necator_rev” denote sets of A. ceylanicum genes with orthology to ES genes in A. caninum, H. contortus, or N. americanus; “ES_Acan.or.Hco.or.L2020” denotes a set of A. ceylanicum genes with orthology to ES genes in any of these three species. “Acan_homology”, “Haemonchus_homology”, and “Necator_homology” denote sets of A. ceylanicum genes with orthology to any genes (ES protein-encoding, or not) in A. caninum, H. contortus, or N. americanus; “Acan.Hco.Nec_homology” denotes a set of A. ceylanicum genes with orthology to any genes in any of these three species. The ES gene sets are either from this study (“all_UMass_ES”) or from Uzoechi et al. [58] (“all_WashU_ES”). Subsets of these ES gene sets that lack orthologies to known ES genes in related parasites are denoted with “non_AcanES_homol_[ES]” (for A. ceylanicum ES genes lacking orthologies to A. caninum ES genes), “non_HcoES_homol_[ES]” (for A. ceylanicum ES genes lacking orthologies to H. contortus ES genes), “non_NecES_homol_[ES]” (for A. ceylanicum ES genes lacking orthologies to N. americanus ES genes), or “non_anyES_homol_[ES]” (for A. ceylanicum ES genes lacking orthologies to ES genes in any of the three other species). Statistical significances of overlaps between orthologous genes and ES genes were computed and described as in S3 Table.(XLSX)
S8 TableA. ceylanicum RNA-seq libraries newly generated in this study.“RNA-seq library ID” provides a short abbreviation for each library; these abbreviations have been used for expression levels (TPMs), mapped read counts (reads), and differential gene expression analyses in later supplementary data tables (S9-S17 Tables). For each library, “SRA accession”, “BioProject accession” and “BioSample accession” provide accession numbers; “Reads”, “Total nt”, and “Read length in nt” give the number of reads, total sequence in nt, and read length; and “RNA-seq description” summarizes biological contents.(XLSX)
S9 TableStatistical analyses of the overlaps between A. ceylanicum genes encoding protein domains (or other traits) and A. ceylanicum genes with either intestine-biased or non-intestine-biased gene expression.Statistical significances of overlaps between domain/trait genes and intestine-biased or non-intestine-biased genes were computed and described as in S3 Table. In addition, pairwise overlaps between intestine- or non-intestine-biased gene sets and other gene categories are provided in the ‘Various’ spreadsheet; since these do not involve multiple comparisons of entire domain sets, only p-values rather than q-values are computed for these overlaps. A. ceylanicum gene categories in “Various” are: “Intest_haemonchus”, genes with orthology to H. contortus genes with intestine-biased expression; “Non-intest_haemonchus”, genes with orthology to H. contortus genes with non-intestine-biased expression; “WashU-intestine-biased”, genes with intestine-biased expression computed from previously published A. ceylanicum RNA-seq data; “WashU-non-intestine-biased”, genes with non-intestine-biased expression computed from previously published A. ceylanicum RNA-seq data (S11 Table) [67,68].(XLSX)
S10 TableThe frequencies with which RNA-seq reads (either from our new RNA-seq libraries listed in S8 Table, or from previously published RNA-seq libraries listed in S11 Table were mapped to A. ceylanicum v2.1 genes by Salmon 1.9.0 [237].(XLSX)
S11 TablePublished RNA-seq data of A. ceylanicum [54,67,68] and H. contortus [185] used in this study.All RNA-seq files were downloaded from the European Nucleotide Archive (ENA). “Species” gives the file’s origin species. “Abbreviation” provides a short abbreviation for each library; these abbreviations have been used for expression levels (TPMs), mapped read counts (reads), and differential gene expression analyses in later supplementary data tables (S12-S17 Tables). “Biological condition (sex, developmental stage, tissue) and replicate number” describes biological content and replication. “Notes” describes any ambiguities that had to be resolved during analysis of the data. “Database” is uniformly ‘ENA’, but is included for completeness in the table. “Accession” and “URL” give SRA accession numbers and ENA Web sources for each file.(XLSX)
S12 TableStatistical analyses of the overlaps between H. contortus genes encoding protein domains (or other traits) and H. contortus genes with either intestine-biased or non-intestine-biased gene expression.Statistical significances of overlaps between domain/trait genes and intestine-biased or non-intestine-biased genes were computed and described as in S3 Table.(XLSX)
S13 TableStatistical analyses of the overlaps between A. ceylanicum genes encoding protein domains (or other traits) and A. ceylanicum genes with positively or negatively immunoregulated intestinal or non-intestinal gene expression.Statistical significances of overlaps between domain/trait genes and positively immunoregulated intestinal genes were computed and described as in S3 Table. In addition, pairwise overlaps between immunoregulated gene sets and other categories (such as “Male-biased”) are provided in the “Various” spreadsheet; since these do not involve multiple comparisons of entire domain sets, only p-values rather than q-values are computed for these overlaps.(XLSX)
S14 TableStatistical analyses of the overlaps between A. ceylanicum genes encoding protein domains (or other traits) and A. ceylanicum genes encoding ES proteins that also have positively immunoregulated gene expression in 19-day intestines.Statistical significances of overlaps between domain/trait genes and positively immunoregulated ES genes were computed and described as in S3 Table. In addition, pairwise overlaps between the positively immunoregulated intestinal ES gene set and other categories (such as “Male-biased”) are provided in the “Various” spreadsheet; since these do not involve multiple comparisons of entire domain sets, only p-values rather than q-values are computed for these overlaps.(XLSX)
S15 TableSoftware used in this study.“Software” provides the name of each computer program or suite of computer programs; “Purpose” describes why this software was used here; “Main web site (URL) or code location” gives the primary Web site for this software; “Bioconda source (if used)” gives the bioconda web site for programs that were installed as bioconda environments; “Online documentation” gives the web site for detailed online manuals, if a given program has one. Publications and arguments for each program are cited and described in Methods.(XLSX)
S16 TablePublished genomic data used in this study.Genome sequences, coding sequences, proteomes, and gene annotations of relevant nematodes were downloaded from WormBase [69] or WormBase ParaSite [194,206] and used for transcriptomic or proteomic analyses; they include Coghlan gene predictions and their corresponding A. ceylanicum genome assembly, along with UniProt [234] and RefSeq [235] proteome databases from highly GO-annotated model organisms. Four separate spreadsheets are given for data files of genomes, proteomes, CDS DNA sets, and gene annotations in GFF format. For each data file, “Species” gives the biological species described by the file, “Comments” describes the particular use(s) to which the data file was put in this study, and “URL” gives the Web source of the file.(XLSX)
S17 TableLC-MS/MS analyses of two A. ceylanicum ES protein sets, E20201108-05 and E20201108-07.Data columns were generated by Proteome Discoverer 2.4; other details of the analysis are given in Methods.(XLSX)
S1 FigGene expression in *A. ceylanicum.*Gene activity is shown for all 33,190 A. ceylanicum genes, with biological replicates on the x-axis and individual genes on the y-axis. Expression levels are in TPM (log10). In S1A Fig, replicates are ordered as in Fig 4; in S1B Fig, replicates are ordered by their similarity of expression (and a dendrogram for their similiaries is shown along the top x-axis). Genes sharing similar patterns of expression are split into 10 clusters. Biological replicates are: young 12-day hookworm adults from normal hosts (YA_12dpi_noDEX); young 12-day adults from immunosuppressed hosts (YA_12dpi_DEX); dissected 19-day intestines from mature hookworm adults in normal hosts (intestines_19dpi_noDEX); dissected 19-day intestines in immunosuppressed hosts (intestines_19dpi_DEX); dissected 19-day non-intestinal tissues in normal hosts (non.intestines_19dpi_noDEX); dissected 19-day non-intestinal tissues in immunosuppressed hosts (non.intestines_19dpi_DEX); previously published adult male intestine (WashU_male_intestine); published adult males (WashU_adult_male); and published adult females (WashU_adult_fem). These heatmaps include two biological replicates (19-day nonDEX intestine replicate 2 and 19-day DEX intestine replicate 2) that clustered anomalously with replicates 1 and 3 of intestinal RNA-seq samples with opposite nonDEX versus DEX status. We interpret this to mean that labeling for the nonDEX/DEX status of these samples was accidentally reversed between RNA harvesting and RNA-seq. We thus omitted these two samples from differential gene expression analysis.(PDF)
S1 FileAn R scrsipt used in batch mode for differential gene expression with edgeR 3.36.0 and R 4.1.3.(TXT)
S2 FileThree R scripts used in batch mode for visualizing differential gene expression in volcano plots with EnhancedVolcano 1.24.0 and R 4.4.3.(TXT)
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Bethony J, Brooker S, Albonico M, Geiger SM, Loukas A, Diemert D, et al. Soil-transmitted helminth infections: ascariasis, trichuriasis, and hookworm. Lancet. 2006;367(9521):1521–32. doi: 10.1016/S 0140-6736(06)68653-4 16679166 · doi ↗ · pubmed ↗
- 2Traub RJ. Ancylostoma ceylanicum, a re-emerging but neglected parasitic zoonosis. Int J Parasitol. 2013;43(12–13):1009–15. doi: 10.1016/j.ijpara.2013.07.006 23968813 · doi ↗ · pubmed ↗
- 3Pullan RL, Smith JL, Jasrasaria R, Brooker SJ. Global numbers of infection and disease burden of soil transmitted helminth infections in 2010. Parasit Vectors. 2014;7:37. doi: 10.1186/1756-3305-7-37 24447578 PMC 3905661 · doi ↗ · pubmed ↗
- 4Beaver PC. Light, long-lasting necator infection in a volunteer. The American Journal of Tropical Medicine and Hygiene. 1988;39(4):369–72. 3189697 10.4269/ajtmh.1988.39.369 · doi ↗ · pubmed ↗
- 5Gems D. Longevity and ageing in parasitic and free-living nematodes. Biogerontology. 2000;1(4):289–307. doi: 10.1023/a:1026546719091 11708211 · doi ↗ · pubmed ↗
- 6Hoagland KE, Schad GA. Necator americanus and Ancylostoma duodenale: life history parameters and epidemiological implications of two sympatric hookworms of humans. Exp Parasitol. 1978;44(1):36–49. doi: 10.1016/0014-4894(78)90078-4 627275 · doi ↗ · pubmed ↗
- 7Ranjit N, Zhan B, Hamilton B, Stenzel D, Lowther J, Pearson M, et al. Proteolytic degradation of hemoglobin in the intestine of the human hookworm Necator americanus. J Infect Dis. 2009;199(6):904–12. doi: 10.1086/597048 19434933 · doi ↗ · pubmed ↗
- 8Maizels RM, Smits HH, Mc Sorley HJ. Modulation of Host Immunity by Helminths: The Expanding Repertoire of Parasite Effector Molecules. Immunity. 2018;49(5):801–18. doi: 10.1016/j.immuni.2018.10.016 30462997 PMC 6269126 · doi ↗ · pubmed ↗