Recruitment of toxin-like proteins with ancestral venom function supports endoparasitic lifestyles of Myxozoa

Annotation and validation of potential toxin-encoding transcripts

Summary statistics and host contamination results for the assembled transcriptomes are given in Table S5. All full-length coding regions were translated, and the predicted amino acid sequences (>30 amino acids in length) were compared to sequences in the UniProt/Swiss-Prot database and our customised toxin dataset (comprised of those in the Tox-Prot database supplemented with putative cnidarian toxin sequences not deposited in UniProt/Swiss-Prot or Tox-Prot databases; see Supplemental Methods). The putative toxins were screened for transcript redundancy and further manually validated by detecting recognised toxin domains and then grouped based on predicted functions (Fig. 2 and Tables S6–S9). A higher number of potential toxins were identified in the C. cruxmelitensis and P. hydriforme datasets in comparison to the myxozoans. Transcripts providing high-scoring BLAST matches to degradative enzymes (e.g. peptidases and lipases in putative toxin protein families) were most diverse in all four datasets, and toxins with predicted haemolytic or coagulation properties were found across all the transcriptomes. The predicted venom transcriptomes of myxozoans contained a notably high proportion of inhibitors of degradative enzymes (especially serine peptidases) (Fig. 2).

The transcriptome of B. plumatellae, obtained from individuals isolated from bryozoan hosts, encoded 124 potential toxins (Fig. 2A and Table S6), the largest group being peptidases (n = 58), including aspartic, serine and metallopeptidases. The metallopeptidases were astacins with and without ShK domains, endoplasmic reticulum aminopeptidases and carboxypeptidases. The next most diverse putative toxin protein family was the serine peptidase inhibitors (n = 23), followed by C-type lectin-like proteins (n = 9), CRiSP - cysteine rich, allergen like proteins (n = 7), and hydrolases including phospholipases (n = 5, with only phospholipase A-like proteins identified), lipases (n = 3), and a phosphatase like protein (n = 1).

The mixed myxosporean transcriptome (Fig. 2B and Table S7) also contained putative toxin peptidases (n = 33), uncharacterised toxin-like proteins (n = 7) and peptidase inhibitors (n = 17, both serine and cysteine inhibitors). We identified cytolytic and coagulation factors including natterins (n = 1) and stonustoxin-like toxins (n = 5), as well as C-type lectin like proteins (n = 4). Fibrogen related proteins (n = 4) and CRiSP—cysteine rich, allergen-like domain containing proteins (n = 3) were found. Toxin-like hydrolases were also identified including lipases (n = 5), phospholipases (n = 5), and a phosphatase (n = 1). One immune factor protein was identified with a homology to alpha-2-macroglobulin. Three transcripts had no recognised domains (see Table S7 for cnidarian toxin homology). It should be noted that as this dataset contained sequences from two different myxozoan species (M. gasterostei and S. elegans), these results represented the diversity of putative toxins shared between both species, rather than the complement of one myxozoan. Although there are some limitations to this interpretation, we surmise that at the level of class (Myxozoa vs other cnidarians) and within the Myxozoa (malacosporeans vs myxosporeans), this dataset provides useful insights into venom toxin diversity.

Transcriptomes of the other two cnidarians contained transcripts for more classical venom component enzymes, with serine proteases and metalloproteinases being the most diverse groups of transcripts. The toxin annotation pipeline identified 302 putative transcripts as potential toxins in the transcriptome of C. cruxmelitensis (Fig. 2C and Table S8). The most diverse putative toxin protein families were the peptidases (n = 121) comprising aspartic, serine and metallopeptidases. Other large groups included fibrogen-related toxins (n = 24, e.g., ficolin or ryncolin-like proteins), phosphatases (n=19), phospholipases including a phospholipase D/DNase like hydrolase (n = 13), and toxin proteins with C-type lectin domains (n = 10). The remaining putative toxins consisted of lipases (n = 8), peptidase inhibitors (n = 5), a “venom component C3” homolog (n = 1) and translationally controlled tumour protein (n = 1). Almost a quarter of the putative toxins (n = 60) had domains found commonly in animal toxins (e.g. SCP/CAP, chitin binding, ShK) but could not be characterised (classed as “Other” in Fig. 2C and see Table S8 for toxin homology). Additionally, 21 putative toxins did not encode any recognised toxin-like domains (7% of the 302).

Peptidases were also the most diverse of the putative toxin protein families annotated from the 214 potential toxin-encoding transcripts in the P. hydriforme transcriptome (Fig. 2D and Table S9). The largest group was the metallopeptidases which included astacins and neprilysin-like proteins. Some metallopeptidase transcripts encoded additional domains in a variety of orders and numbers (e.g. astacins with or without repeats of TSP, MAM and ShK domains). Other putative toxins, such as 5′ nucleotidases (n = 18), lipases (n = 8), both A and B phospholipases (n = 9), oxidases (n = 6), C-type lectin like proteins (n = 8) were also identified. As for C. cruxmelitensis, a large proportion of the putative toxin transcripts encoded recognised toxin domains, but were not characterised further (n = 23, ~13% of the flagged transcripts, see Table S8 for cnidarian toxin matches). Fifteen putative P. hydriforme toxins had no recognised domains. A higher proportion of putative toxin transcripts had predicted N-terminal signal peptides in P. hydriforme (42% (89/214)) than in C. cruxmelitensis (23% (70/302)) (Tables S8 and S9).

Annotation of potential toxin-encoding peptides

Based on analyses of the myxozoan proteomes, only a small number of peptides were detected and identified as putative toxins compared to the number of toxins identified in the myxozoan transcriptomes (Fig. S6). Four of 124 putative toxin transcripts (3%) were translated in B. plumatellae and two of 96 in the mixed myxosporean material (2%) (see Tables S6 and S7). The four putative B. plumatellae toxins were assigned to three toxin protein families (CRiSP, C-type lectin and serine peptidase inhibitor) with potential neurotoxin, dysregulation of blood haemostasis and enzyme inhibitor biological functions (Table S6). The two putative toxins in the mixed myxosporean proteome were linked with similar functions, one being neurotoxin (verrucotoxin) and the other an inhibitor of blood haemostasis (Blarina toxin) (Table S7). These toxins may be shared between both M. gasterostei and S. elegans or unique to only one.

The numbers of potential toxin encoding transcripts translated and detected in the proteomes of C. cruxmelitensis and P. hydriforme were substantially higher than in the myxozoans (Fig. S6). In C. cruxmelitensis 64 of 302 transcripts (21%) were translated and in P. hydriforme 77 of 214 (35%) were translated (Table S10). The predicted toxins of C. cruxmelitensis venom had a broad range of predicted biological activities, the most prominent being cytotoxic (62 of 77 transcripts (81%)). Two putative neurotoxins, a phospholipase A2, and a serine protease inhibitor with homology to actitoxin of Anemonia viridis were detected. Eleven proteins lacked supporting data to ascribe a biological function or were involved in toxin maturation (e.g. glutaminyl-peptide cyclotransferase) (Tables S8 and S10). Peptides that could be annotated with predicted cytotoxic activity were part of either degradative enzymes or inhibitors of haemostasis and included metalloproteases (n = 25, including two proteins similar to snake venom metalloproteases, SVMPs), proteases (n = 9), lipases and phospholipases (n = 9), oxidases (n = 6), hydrolases (n = 5), and C-type lectins (n = 8). Most of the 77 translated putative toxins in the proteome of P. hydriforme (73%) (Tables S9 and S10) were annotated as degradative enzymes (peptidases (n = 33), lipases (n =10), glycoside hydrolase (n = 2) and 5′-nucleotidase (n = 2)) or inhibitors of blood haemostasis (C-type lectins (n = 8)). Peptides with predicted neurotoxin activity were considered minor components of the P. hydriforme venom proteome and included phospholipase A2 (Tables S9 and S10).

Comparative venom profiles based on proteomics data

The presence or absence of translated proteins assigned to toxin families in Anthozoa, Medusozoa, P. hydriforme and myxozoans (Table S1) is summarised in Fig. 3. The Myxozoa had the smallest number of toxin protein families (n = 7), followed by P. hydriforme (n = 15), Anthozoa (n = 29) and Medusozoa (n = 29). Four toxin protein families were common to all groups (C-type lectins, serine peptidases, serine peptidase inhibitors, and CRiSP proteins). Six toxin protein families were unique to Medusozoa or Anthozoa, and an additional eight toxin protein families were shared only between these two lineages. No toxin protein families were unique to Myxozoa. Myxozoans shared six of their seven toxin protein families with both Anthozoa and Medusozoa, and one uniquely with Medusozoa. Polypodium hydriforme shared eight toxin protein families with Medusozoa and 10 jointly with Medusozoa and Anthozoa. Aspartyl peptidase was unique to P. hydriforme, having never been described previously in the predicted venom of any cnidarian (Fig. 3).

Evolutionary analysis

We used a molecular phylogenetic framework to clarify evolutionary trends and investigate the distribution of toxin-like genes found in Endocnidozoa and other animals. Phylogenetic hypotheses of the relationships between amino acid sequences of the putative venom Kunitz-type and phospholipase A2 were inferred in this study from a wide taxonomical range of venomous and non-venomous animals possessing these gene families (Tables S2 and S3). These two families were chosen because they have been frequently studied in the cnidarian literature and other venomous lineages, allowing reasonable comparative standards. Phylogenetic trees constructed from multiple alignments of amino acid sequences (Figs. S4 and S5) using maximum likelihood (venom Kunitz-type is shown in Fig. 4 and phospholipase A2 is shown in Fig. 5) and Bayesian inferences (Figs. S7 and S8) for both protein families have similar topology.

Phylogenetic analysis of the venom Kunitz-type protease inhibitors revealed an evolutionary pattern only partially congruent with the expected taxonomy, in which more inclusive groups (collapsed in Fig. 4) have low bootstrap support (≤60%). Cnidarian sequences are split in several less inclusive groups or as terminals along the tree (Fig. 4 and Table S11). Cnidarian groups were not resolved as monophyletic (the exception being Staurozoa whose outcome depended on resolution of the polytomy). Three endocnidozoan groups were supported (2 for Buddenbrockia plumatellae and 1 for mixed myxozoans). At least one sequence of each group of Endocnidozoa (i.e. B. plumatellae, mixed myxozoans, and Polypodium hydriforme) grouped with a venomous taxon (scorpions, spiders and conid snails), and, they never grouped together with other cnidarians. The substitution rates of the amino acid sequences of the venom Kunitz-type are significantly different when comparing Anthozoa and Endocnidozoa (P < 0.05; P = 0.023), with endocnidozoan lineages evolving slightly faster (cf. RRTree version 1.1; (Robinson-Rechavi & Huchon, 2000)) than those of free-living cnidarians.

The high bootstrap values in the phospholipase A2 tree (Fig. 5) support three distinct groups. Sequences in Group I (in red) and Group II (in green) derived from venomous snakes and non-venomous mammals. The third group (Group III; in blue) contains a basal sub-group composed of anthozoans which is sister to a sub-group comprising medusozoans (all four classes) and another sub-group with sequences of medusozoans, anthozoans, venomous non-cnidarian taxa (three scorpions, one honey bee, one lizard) and a sequence of B. plumatellae.

Localisation of putative toxins in B. plumatellae

Two out of the three polyclonal antibodies raised against custom synthesised putative toxins identified from both transcriptomic and proteomic data for the myxoworm B. plumatellae isolated from bryozoan hosts were localised in early presporogonic stages (worms lacking spore development) (Figs. 6 and 7). Fluorescence corresponding to antibodies raised against a C-type lectin (Fig. 6) and to a serine peptidase inhibitor (Figs. 7C, 7D) was observed in the cytoplasm of cells internal to the muscle layer. Sporogony in B. plumatellae occurs internal to the muscle sheet, filling the inner cavity with malacospores (see Fig. 7A). Cells in this region are therefore assumed to be presporogonic stages. Both C-type lectin- and serine peptidase inhibitor expression was also observed in early sporogonic cells (Figs. 6C and 6D; Figs. 7C and 7D). Additionally, the serine peptidase inhibitor antibody was expressed in the cytoplasm of outer epithelial cells (Figs. 7B) and possibly on the external surface of malacospores (Figs. 7E and 7F). No signal was observed for the putative B. plumatellae CRiSP allergen in any of the developmental stages examined.

Conservation, co-option and divergence of toxins

We present evidence that toxin genes have been conserved and repurposed in endoparasitic cnidarians, in support of H1. In particular, multiple bioinformatic analyses of transcriptomic and proteomic data infer that endocnidozoans utilise proteins with sequence homology to toxins found in venoms of free-living cnidarians and other animal taxa. In addition, localisation studies provide evidence for expression of selected toxins in cell populations involved in spore development and sub-epithelial cells but not within nematocyst-containing spores. The latter observation suggests a function other than venom delivery but further study is required to discount alternate explanations such as lack of detection (e.g. ultrastructural investigations to localise toxins in nematocysts or in sporoplasmosomes in spore cells), stage of maturity (e.g. toxins not yet transferred), or lack of dye penetration. Nevertheless, demonstration of toxin expression in itself illustrates a repurposing of venom activity because it is not associated with prey capture or defence as in free-living cnidarians. Notably, P. hydriforme achieves transmission using a specialised gonadophore that has been observed to attach to juvenile (pre-larvae) sterlet with nematocysts (Raikova, 1994) but whether toxins are involved is unknown.

Molecular phylogenetic analyses of the two expressed protein families show that the Endocnidozoa sequences may be closer to sequences from other cnidarians (like for PLA2, Fig. 5), although the molecular phylogenetic analyses of the Kunitz-type amino acid sequences have not elucidated the evolutionary relationships between both groups (Fig. 4). The resolution of the topology recovered with Kunitz-type is low and mainly supports deep nodes. Low bootstrap values can indicate conflicting or little signal in the data set. The latter may arise when sequences are very conserved, resulting in a data set with few informative sites and lack of information—a problem that can be exacerbated when analysing a limited and small functional domain such as Kunitz-type. Such sequence homogeneity can occur through processes such as natural selection, randomness or by intrinsic biases in the production of the molecular variability (Storz, 2016).

Interestingly, however, some endocnidozoan Kunitz-type sequences of B. plumatellae, mixed myxozoans and P. hydriforme group with sequences of venomous animals (other than cnidarians) providing evidence for divergence of venoms from ancestral free-living cnidarians (addressing H2) and potential convergent evolution of venoms in diverse animal groups. This contrasts with other myxozoan inhibitors (non-toxic) that group by taxonomy rather than physiological role (Eszterbauer et al., 2020). Fuzzy relationships with low bootstrap values and poor resolution were corroborated in a much broader analysis encompassing 356 sequences presenting the BPTI/Kunitz inhibitor mature domain obtained from UniProtKB/Swiss-Prot and including several families of Kunitz-type proteins (Figs. S9 and S11). The various phylogenetic analyses of the Kunitz-type toxin family suggested sequences were highly conserved which may account for lineage-specific gene duplications in many groups, including Endocnidozoa. Duplications may enhance the number of paralogs that then diversify in association with different biological functions, potentially further supporting H2.

Endocnidozoan toxin diversity

The transcriptome for B. plumatellae stages developing in invertebrate hosts revealed an array of toxin-encoding transcripts similar to venom toxins of free-living cnidarians and a notably high percentage predicted to encode for serine protease inhibitors. However, relatively fewer toxin-encoding transcripts were demonstrated in myxozoans (124 in B. plumatellae; 96 shared between the mixed myxosporeans) than in C. cruxmelitensis (302, Fig. 2) and other free-living cnidarians. For example, some 170, 218 and 508 potential toxin transcripts were identified in transcriptomes of the cubozoan Chironex fleckeri, the scyphozoan Stomolophus meleagris, and the anthozoan Stichodactyla haddoni, respectively, on the basis of homology to known toxins in publicly available sequence databases (Brinkman et al., 2015; Madio, Undheim & King, 2017; Li et al., 2014). The detection of fewer toxin genes in myxozoans supports H3 and is concordant with reduced genomes sizes of myxozoans (Chang et al., 2015). In this case, loss may relate to relinquishing a predatory lifestyle. No myxozoan specific putative venom toxins were identified in the mixed myxosporean transcriptome or proteome datasets.

The predicted venom proteome of P. hydriforme contained a comparable number of putative toxins to those in the C. cruxmelitensis venom proteome (77 and 64, respectively). These results and the lack of detection of novel toxins in B. plumatellae suggest that toxin diversity is not linked with parasitism per se. In addition, they are consistent with evidence that life cycle complexity is not associated with an increased number of genes in cnidarians (Gold et al., 2019) (P. hydriforme has both free-living and parasitic stages while C. cruxmelitensis has a free-living, benthic stage only).

Toxin demonstration, translation and function

Adopting a stringent approach to annotate putative toxins in transcriptomes and proteomes is crucial because toxin encoding gene families derive from housekeeping gene families. Expression of non-toxin homologues makes it notoriously difficult to identify toxins in homology-based annotations (Li et al., 2014; Hargreaves et al., 2014; Macrander, Broe & Daly, 2016). One approach to confirm toxin identity is to assess whether high expression levels are associated with venom glands (Smith & Undheim, 2018). This approach is precluded for Cnidaria which package toxins into intracellular organelles (nematocysts). Consequently, methods to distinguish transcripts encoding toxin from non-toxin homologs using sequence and domain homology searching methods are now standard in cnidarian venom research and were included here (see (Macrander et al., 2018) and references therein). In addition, we undertook molecular phylogenetic analyses of selected putative venom toxins analysing proteins in both venomous and non-venomous animals. The clustering of mainly Endocnidozoa sequences in both trees with homologs from other venomous animals provided further evidence that the sequences are toxins. These combined approaches greatly strengthen the evidence for the presence of genes that encode for venom toxins. Further substantiation of this scenario could include support for involvement in host interactions (e.g. localisation to excretory vesicles (sporoplasmosomes)) or demonstration that housekeeping homologs in venom-free tissues are distinct from putative toxins in molecular phylogenetic analyses.

There were surprisingly few matches between proteomic and transcriptomic sequences in the myxozoan data compared to their free-living counterparts (Fig. S6). Although 124 transcripts were predicted to encode for toxins in B. plumatellae, only four were translated (3%; Table S6). Similarly, only 2 of 96 mixed myxosporean transcripts (2%; Table S7) matched to the proteome data highlighting the importance of using proteomics to validate identification of putative toxins from transcriptome sequences, whilst also minimizing false positive identification of non-toxin homologs identified in transcriptomes (Madio, Undheim & King, 2017). Myxozoan proteomes revealed toxins in families with variable functions. These included procoagulant and anticoagulant toxins (C-type lectins in B. plumatellae; function demonstrated in snakes (Morita, 2004)) and neurotoxins such as verrucotoxin and blarina-like toxin from the mixed myxosporeans (shown to have vasodilation and paralysis properties in fish and shrews (Kita et al., 2004; Yazawa et al., 2007)). Additionally, cysteine-rich secretory proteins (CRiSPs) and serine peptidase inhibitors identified in B. plumatellae are associated with neurological and muscular dysfunction as ion channel blockers (Harvey, 2001; Yang et al., 2014) and with immune system modulation (Falcón et al., 2014; Ranasinghe et al., 2015).

Our localisation studies indicate translation of putative toxins in B. plumatellae. We found both the C-type lectin and serine peptidase inhibitor in presporogonic cells. Some epithelial cells forming the outer wall of a myxoworm also exhibited expression of serine peptidase inhibitor suggesting the possibility of toxin redeployment, for example in the breakdown of the muscle tissues (Canning et al., 2002) enabling spore release. There was also evidence for serine peptidase inhibitor on the outside of mature malacospores. These collective observations suggest that venom toxins are uniquely deployed in these endoparasites and further localisation studies in different life stages are now warranted to investigate possible repurposing of toxin function.

In the case of P. hydriforme there was a high proportion of matches between the proteome and transcriptome datasets (Table S9). Degradative enzymes (peptidases, lipases, hydrolases) made up the largest complement of putative toxins, representing 62% (48/77) of the toxin repertoire. These enzymes may assist the free-living stolon stage of P. hydriforme to digest remaining yolk reserves that fuel adoption of benthic life or they may be used for prey capture once yolk reserves are depleted and a mouth has formed (Raikova, 1994).

Expanded understanding of cnidarian toxins and toxin families

Our three endocnidozoan venom proteomes expand the current picture of cnidarian venom composition while our collective evidence suggests that myxozoans produce complex and varied toxins that could have functioned as venoms in free-living ancestors. We observed one such toxin in cell populations of the outer body wall and the same toxin and another in cells involved in spore development. The former observation suggests repurposing of function from prey capture/defence via venom-delivering nematocysts to another function such as tissue degradation that facilitates release of spores from myxoworms (Canning et al., 2002). Whether toxins produced during spore development are eventually stored in nematocysts within spores requires further investigation.

Proteomic data suggest that 23 toxin protein families are shared amongst cnidarians (Fig. 3) and patterns of sharing between Endocnidozoa and other cnidarians imply an early diversity and long evolutionary history of a common venom complement (Jaimes-Becerra et al., 2019; Doonan et al., 2019). Notably, with the exception of ficolin lectin, all of the putative toxins detected in the myxozoan transcriptomes and proteomes have been reported variously in anthozoans (Fig. 3, Tables S1 and S6–S9). Our data thus provide evidence for the conservation of many core toxin families across free-living and endoparasitic cnidarians and raise questions about the early evolution of the venom trait.