Next Article in Journal
Individual and Combined Effects of Aflatoxin B1 and Sterigmatocystin on Lipid Peroxidation and Glutathione Redox System of Common Carp Liver
Previous Article in Journal
Cytotoxic and Genotoxic Effects of Cyanobacterial and Algal Extracts—Microcystin and Retinoic Acid Content
Article

Diet Diversity in Carnivorous Terebrid Snails Is Tied to the Presence and Absence of a Venom Gland

by 1,2,3,4,†, 1,5,6,†, 1,7 and 1,2,3,*
1
Department of Chemistry, Hunter College Belfer Research Center, City University of New York, New York, NY 10021, USA
2
Graduate Programs in Biology, Biochemistry, Chemistry, Graduate Center, City University of New York, New York, NY 10016, USA
3
Division of Invertebrate Zoology, The American Museum of Natural History, New York, NY 10024, USA
4
Department of Biology, Hofstra University, Hempstead, NY 11549, USA
5
Department of Biology and Biotechnologies “Charles Darwin”, Sapienza University of Rome, I-00185 Rome, Italy
6
Department of Biology and Evolution of Marine Organisms (BEOM), Stazione Zoologica Anton Dohrn, I-00198 Rome, Italy
7
Department of Ecology, Evolution, and Marine Biology, University of California Santa Barbara, Santa Barbara, CA 93106, USA
*
Author to whom correspondence should be addressed.
These authors contributed equally.
Received: 11 December 2020 / Revised: 21 January 2021 / Accepted: 26 January 2021 / Published: 2 February 2021
(This article belongs to the Section Animal Venoms)

Abstract

Predator-prey interactions are thought to play a driving role in animal evolution, especially for groups that have developed venom as their predatory strategy. However, how the diet of venomous animals influences the composition of venom arsenals remains uncertain. Two prevailing hypotheses to explain the relationship between diet and venom composition focus on prey preference and the types of compounds in venom, and a positive correlation between dietary breadth and the number of compounds in venom. Here, we examined venom complexity, phylogenetic relationship, collection depth, and biogeography of the Terebridae (auger snails) to determine if repeated innovations in terebrid foregut anatomy and venom composition correspond to diet variation. We performed the first molecular study of the diet of terebrid marine snails by metabarcoding the gut content of 71 terebrid specimens from 17 species. Our results suggest that the presence or absence of a venom gland is strongly correlated with dietary breadth. Specifically, terebrid species without a venom gland displayed greater diversity in their diet. Additionally, we propose a revision of the definition of venom complexity in conoidean snails to more accurately capture the breadth of ecological influences. These findings suggest that prey diet is an important factor in terebrid venom evolution and diversification and further investigations of other understudied organisms, like terebrids, are needed to develop robust hypotheses in this area.
Keywords: diet diversity; venom; Terebridae diet diversity; venom; Terebridae

1. Introduction

Venomous animals use their toxin arsenal for several ecological activities including predation, defense, space competition, and intraspecific communication [1]. In the predator–prey arms-race, diet is widely accepted as a major force in shaping venom evolution across different taxa such as scorpions, snakes, and spiders [2,3,4]. However, there have been few multi-species studies that specifically look at the role of diet on venom composition and complexity [5,6]. Recently, dietary breadth of venomous snails in the family Conidae was positively linked to venom complexity (number of compounds in venom), but not venom composition (type of compounds found in venom) [7]. Such results give credence to a possible association between venom evolution and prey acquisition throughout the animal kingdom [8,9]. Venom, a cocktail of toxic compounds that includes peptides, proteins, and small molecules, is produced to elicit a reaction with targeted organisms. Since venom is fundamentally a molecular trait, the correlation between venom composition and prey phenotype can be used to trace adaptive evolution and ecological processes such as predation and diet [1,10].
There are primarily two current hypotheses to explain the relationship between a venomous predator’s diet and the nature of its venom: (1) prey preference can determine venom components (types of compounds in venom) [4], and (2) dietary breadth and venom complexity (number of compounds in venom) are positively correlated [11,12]. Studies showing that venoms from different species are most effective on their preferred prey help support the first hypothesis that prey preference can drive venom composition [13,14]. Similarly, to support the second hypothesis (the number of compounds in venom), studies have shown that increased prey breadth is correlated with venom complexity. Examples supporting each hypothesis include findings by Phuong et al. (2016) that cone snails with a more generalized diet tend to have more complex venoms [7] and Li et al. (2005) who found a 50- to 100-fold decrease in toxicity of venom in a group of sea snakes that have shifted to eating only fish eggs [7,11]. Knowledge of selective forces, such as diet, that can drive venom diversification is crucial to understanding patterns in venom evolution.
A working theory in venom evolution is that the acquisition of a venom apparatus and the capacity to synthesize venom and deliver venom to prey has led to diversification within clades of venomous organisms. The acquired novelty of venom then presents opportunities to feed on more or different types of prey, thereby opening new niches and facilitating speciation [15,16]. In contrast, we recently found no correlation between the presence of a venom gland and increased diversification rates within the Terebridae (auger snails) [17]. Terebrids are an understudied family of carnivorous marine snails whose structural variations in the foregut anatomy denote a variety of prey capture techniques. In this study, we address the hypothesis that repeated innovations in terebrid foregut anatomy and in their venom composition may correspond to variation in the diet of terebrids. [18,19]. Similar to cone snails, certain species of terebrids produce a complex mixture of peptide and protein toxins in their venom to subdue their annelid prey [20,21,22,23]. Several members of the Terebridae deliver their toxin payload using a venom apparatus comprised of a venom gland, proboscis, and radula sac, similar to the venom apparatuses found in cone snails. Although the ancestral state for terebrid anatomy included a venom apparatus, many lineages of terebrids have since lost a structured venom apparatus and we recently reported twelve distinct terebrid foregut anatomies, only five of which have a traditional venom gland [17]. Those species that lost the structures of the venom apparatus may still possess other accessory organs, such as salivary glands and proboscis feeding structures [17]. Our previous work also found that each foregut type could be found across multiple clades of the terebrid phylogeny [17].
Here, we examine the correlation between terebrid foregut anatomies (that imply the presence or absence of a venom apparatus), venom complexity, and prey diversity, by metabarcoding the gut contents of 71 specimens sampled across 17 terebrid species (Figure 1, Supplemental Table S1). The terebrid species selected span eight terebrid genera, belong to five phylogenetic clades (clades B, C, D, E, and F) out of the total six found in the family, and are representative of at least four of the twelve known terebrid foregut anatomies [17]. We used gut and venom gland transcriptomes from each of the 71 terebrid specimens to identify the annelid worm meal, the corresponding venom composition, and to examine the role of diet in the diversity of terebrid venom apparatus and venom composition.

2. Results

2.1. Terebrid Gut Content Reveals Diverse Annelid Prey

We sequenced the gut content of 17 snail species, belonging to 5 terebrid clades (clades B, C, D, E, and F). The bioinformatic pipeline firstly isolated a total of 254 prey reads (Table S2). After merging data from species replicates, we obtained 217 unique record indicating the feeding of one snail species on one prey genus (Tables S1 and S2). The range of worm prey genera detected in each terebrid specimen was between 1 and 7 (M = 3.06, SD = 1.37). Worm meals were identified in terebrid lineages with and without a venom apparatus. The range of worm prey genera detected for each terebrid species was between 2 and 11 (M = 6.18, SD = 2.58) (Figure 2).
Overall, eighteen different worm genera, belonging to nine families and two phyla (Anellida and Sipuncula) were identified as terebrid prey. The most represented annelid prey was from the genus Scolelepis Blainville, 1828 (Spionidae: Polychaeta), which was found in 16 species out of 17 terebrid species and in 89% of the analyzed specimens (n = 71). Three other annelid genera, Pseudopolydora Czerniavsky, 1881 (Spionidae: Polychaeta), Spio Fabricius, 1785 (Spionidae: Polychaeta) and Neanthes Kinberg, 1865 (Nereididae: Polychaeta) were also largely represented in our dataset and were found, respectively, in 13–15 terebrid species and 30–48% of our specimens. DNA of each of these four prey genera were detected in terebrid species collected in all macro sampling areas (Florida, UAE, and Papua New Guinea), in terebrids with and without a venom apparatus, and belonging to five phylogenetic clades investigated. In general, the Annelid genera found as terebrid prey belonged to two classes: Errantia (7 genera) and Sedentaria (10 genera). In all specimens but one, we found at least one Sedentaria prey and all other samples had between 1 and 5 Sedentaria prey genera, while in only half of the dataset (n = 36) we detected at least one Errantia prey, with the highest number being 4 (Figure 2). Additionally, a live Terebra guttata sample (KVG_256), collected in Papua New Guinea, was found with a worm in its mouth (Figure 1). We sequenced the 16S rDNA gene of the worm (GenBank accession number MW013507) and BLASTed it against the non-redundant protein sequences database. This approach identified the annelid species as Scolelepis eltaninae with a reported e-value of 7 × 10−94. S. eltaninae belongs to the Spionidae family, which was highly prevalent in all terebrid species analyzed in this study. Su

2.2. Strong Correlation between Presence and Absence of a Venom Apparatus and Terebrid Diet Diversity

The Kruskal Wallis non-parametric test was performed on all 71 specimens and the result indicated a significant difference of diet breadth between terebrid specimens with and without venom apparatus (p < 0.01). Specifically, terebrid specimens with a venom gland present had a more restricted diet (12 annelid genera recorded) than terebrid species without a venom gland (16 annelid genera recorded) (Figure 3A). A randomization test of all 71 specimens was conducted to further corroborate the results of the Kruskal Wallis non-parametric test. The difference between the Shannon diversity indices observed for the samples with the presence (1.91) and absence of a venom apparatus (2.26) was calculated and found to be 0.352. Observations were then randomized by column 100 times and the differences in Shannon indices were recalculated. Our data produced the greatest difference in Shannon index (p < 0.01), supporting the finding that terebrid species without a venom gland had a greater diversity of annelid diet (Figure 3B).

2.3. No Correlation between Terebrid Venom Complexity and Diet Diversity

The transcriptomes from six terebrid specimens (Hastula hectica, Hastula matheroniana, Terebra subulata, Terebra guttata, Terebra argus, and Terebra straminea) were compared with the results of their gut content eDNA to investigate correlations between terebrid venom complexity and diet diversity, where venom complexity was defined by the number of total putative terebrid toxins and the number of cysteine frameworks identified. The number of mature terebrid venom peptide toxins (teretoxins) was quite variable, ranging from 41 to 128 in our analyzed species (Figure 4A).
A small variation in the number of cysteine frameworks, ranging from 8 to 10 cysteines, was instead detected among the examined terebrid species (Figure 4B). Hastula hectica showed the lowest venom complexity while Hastula matheroniana showed the highest. Plots of the average species Shannon index against the number of teretoxins and cysteine frameworks did not produce any quantitative statistical results, indicating there is no direct correlation between terebrid venom complexity and diet (Figure 4).
However, there are qualitative results worth noting. For example, Hastula hectica has a venom composition different from the other five individuals analyzed in this study (Figure 5).
Specifically, H. hectica has a higher proportion of frameworks I and XI putative teretoxins and surprisingly has no framework XXII putative teretoxins, which on average account for approximately 7% of all teretoxins found in all specimens tested. Interestingly, of the six individuals analyzed, H. hectica was the only individual that consumed Neanthes and did not consume Scolelepis, which was found as prey in all other terebrid species examined (Figure 5). Likewise, T. subulata, T. argus, and T. straminea all consume the same genera of annelids (Scolelepis and Pseudopolydora) but their venom composition is not the same. A similar pattern is shown with H. matheroniana and T. guttata, further suggesting venom complexity may not be associated with diet preferences.

2.4. No Correlation between Shell Length and Diet Diversity

To investigate other factors of terebrid diversification that might not be tied to venom, we investigated abiotic characters pertaining to shell length, geographic location, and depth using all 71 specimens. No correlation emerged between shell length and number of prey genera for each specimen (Figure S1B). The Jaccard heat map identified similarity in diet between Myurella undulata and Terebra dislocata. However, Shannon’s diversity index (H’) distributions showed nonsignificant differences between phylogenetic clades, macro collecting areas, and shell lengths (Figure S1A). In addition, no significant differences were also detected when comparing H’ distributions and maximum, minimum, average and range of collection depth. These findings may suggest abiotic pressures may have a limited role in terebrid annelid diet diversity.

3. Discussion

How the diet of venomous animals influences the composition of venom arsenals is an ongoing debate [7,17,24]. Elucidation of the selective forces that drive diversification of venomous organisms, requires examination of variations in mechanisms of delivery (beaks, fangs, harpoons, spines, etc.), disparities in venom composition (proteins, peptides, small molecules, etc.), and prey preferences. The degree of anatomical variation found in the foregut of terebrids, including the presence or absence of a comprehensive venom delivery apparatus, was previously thought to be correlated with venom composition and annelid prey preference [18,25,26]. We examined this hypothesis using a next generation sequencing approach, molecularly identifying (for the first time) prey genera in the gut contents of terebrids. We identified a strong correlation between the presence and absence of a venom apparatus and the diversity of the terebrid diet.
In this study, we detected a total of 18 annelid prey genera being consumed by terebrids, of which Scolelepis and Spio (both in the Spionidae family) were the most prevalent (Figure 1 and Figure 2). Using these data, we found terebrid species with a venom apparatus consume fewer annelid prey types than those without a venom apparatus. Our results were further corroborated with a randomization test and are not an artifact of the sampling group. Our results correspond with the current literature that evolution of venom production can cause prey specificity [7,9,27,28]. In a recent study using snakes, it was found that venom diversity was a result of adaptations to specific diets and the complexity of venom in snakes is due to the fact that venoms are under strong natural selection pressures to combat selection pressures in prey molecular targets [4,27,29]. As with snakes, the qualitative venom variation observed in terebrids with a venom gland (specifically Hastula hectica) might be explained by an evolutionary arms race between a terebrid predator and a specific annelid prey, however, additional taxa are needed to corroborate this finding for terebrids.
While the diversity of foregut anatomies found in the entire Conoidean superfamily is present in just the Terebridae family, contrary to what was found for cone snails, we did not find a statistical correlation between venom complexity and diet diversity [7,30]. The lack of statistical significance may reflect: (1) the small sample size (only 6 of the 17 terebrid species were used in the venom variation analysis) and (2) a need to redefine venom complexity. In Phuong et al.’s paper (2016), it was determined that total number of mature toxins was the best indicator of venom complexity, while the number of gene superfamilies and the number of cysteine frameworks were less ideal indicators [7]. Here we used both the total number of mature teretoxins and the diversity of cysteine frameworks to determine venom complexity but did not find a statistical correlation. Although there was no statistical correlation between the genera of prey consumed and the distribution of cysteine frameworks in the venom, there is a qualitative relationship wherein Hastula hectica had a different suite of putative teretoxins than the other terebrids analyzed. H. hectica was the only terebrid found to consume Neanthes and not Scolelepis. From our analyses we cannot confirm if the diet preference in H. hectica can be explained by the different putative teretoxins expressed in its venom gland. However, contrary to other terebrids examined, H. hectica are found in the intertidal zone in the surf. Their restricted habit could be why they are the only ones that feed on Neanthes and not Scolelepis which are also found in the surf. H. hectica results highlight how important it is to examine teretoxins and diet within the context of a single individual. The tremendous amount of variation in terebrid foregut anatomy can possibly be explained by the niche variation hypothesis. The niche variation hypothesis indicates “[among-species] populations with wider niches are more variable than populations with narrower niches” [31]. The terebrids, which are globally distributed and found at bathymetric sea levels ranging from intertidal to up to the pelagic zones, have broadened their niches and their diet so that they can overlap in geographical location and possibly even niche, but still do not have to over-compete for food [26]. Notably, we recently determined environmental pressures may have more of an influence on terebrid diversification than the presence of venom [17].
Overall, we determined that terebrid diet breadth can be correlated with the presence or absence of the venom gland. Specifically, those terebrid species used in this study with a venom gland showed a more restricted diet than those without it. On the other hand, our data suggest that other factors, such as terebrid species classification, shell size, collecting depth, biogeography, and venom complexity, do not significantly influence diet diversity. Additionally, if we consider also the qualitative data we collected, these suggest that, apart from H. hectica, terebrids overlap in the prey genera that they consume. Our findings highlight several possibilities about terebrid venom evolution, such as: (1) terebrid species that lack a venom gland may still produce venom via a different tissue, such as salivary glands, allowing them to capture and consume the same genera of annelid prey. It has been reported that the salivary glands of members of the family Conidae (e.g., Conus pulicarius), phylogenetically close to Terebridae, can produce venom toxins, as well as those of more distant taxa such as leeches, ticks, and others [32,33,34,35,36,37]. Indeed, several terebrid species used in our study lack a venom gland but have salivary glands. (2) Terebrids without any components of the venom apparatus (venom gland, proboscis, radular sac) may either release venom directly into the water, as has been document from cone snails, or may forego using venom to subdue their prey. If this is true, it could explain why the traditional venom gland has been lost on eight occasions in the Terebridae [19]. Taken together, our findings and those of Phoung and colleagues suggest the current method of classifying cono- or tere-toxins is lacking in details and proposes that venom composition in conoidean snails should be redefined to more accurately reflect the complexity. A holistic model in which we examine, both number of toxins, cysteine framework, gene superfamilies, and prey composition may provide more valid analyses of venom arsenal classifications. Finally, this work illuminates the importance of examining neglected venom organisms in constructing robust theories about venom evolution and species diversification.

4. Materials and Methods

4.1. DNA Extraction and Amplification for Next Generation Sequencing

Seventy-one terebrid species (n = 71) collected in Kavieng, Papua New Guinea (n = 60), Tampa, Florida (n = 4), and Fujairah, United Arab Emirates (n = 7) were chosen from five phylogenetic clades to represent the diverse foreguts found in the Terebridae family. DNA was extracted using the following phenol chloroform extraction protocol. Dissected guts were lysed overnight in lysis buffer and proteinase K at 56 °C. After centrifuging the samples post incubation, 200 µL phenol was added. Samples were centrifuged at 13,000 rpm and the supernatant was added to 100 µL phenol and 100 µL chloroform. Following centrifugation, the supernatant was added to 200 µL chloroform and centrifuged. The supernatant was added to 20 µL sodium acetate (3M, pH 5.2) and 400 µL 95% ethanol. After 2 h resting at −20 °C, samples were centrifuged and the resulting pellet was washed in 500 µL 70% ethanol. The pellet was collected and suspended in 50 µL Tris EDTA buffer (pH 8) and stored at −20 °C. DNA was quantified usinga Qubit fluorometer (Invitrogen, Carlsbad, CA, USA).

4.2. PCR Primers and Amplification

We established a set of PCR primers and a PCR protocol that optimized the amplification of the 16S rDNA gene in marine worms. The PCR primers should ideally find a balance between the short, degraded prey DNA found in the guts and a viable sequence length necessary for identifying the prey. We aimed to optimize primers that would amplify a wide range of marine worm DNA, the known prey of terebrids, without amplifying the terebrid DNA. For this study, we amplified a ~350 bp fragment of the 16S rDNA gene using a pair of primers designed to amplify worms DNA and to exclude that from the mollusc: 16SANNF3 (GTATCCTGACCGTGCWAAGGTAGC) designed by 37, and 16Spr1_NP (CCTAAGCCAACATCGAGGTGC) which was modified from the primer 16Spr1 by 37 to not amplify turrid 16S rDNA as well (N. Puillandre, pers. com.). The efficiency of primers was tested by amplifying a mix of DNA containing terebrid and marine worm DNA in the ratios 100:1 and 1000:1, respectively. PCR products were confirmed on a 1.3% agarose gel stained with GelRed.
PCR was carried out in 25 µL volumes with a final concentration of 1× of the supplied PCR buffer, Platinum Hot Start PCR Master Mix 2X (Invitrogen, Carlsbad, CA, USA) and 0.2 µM of each forward and reverse primer. For each reaction, 1 µL of the extracted DNA was added. The thermocycler protocol was 94 °C, followed by 40 cycles of 94 °C for 20 s, 60 °C for 20 s, and 72 °C for a minute. The extension period lasted for 5 min at 72 °C and the PCR products were cooled at 17 °C and stored at −20 °C following PCR product confirmation. Confirmation of PCR products was based on the presence of a band in a 1.3% agarose gel stained with GelRed.

4.3. Library Preparation and Sequencing

Libraries were prepared using the Nextera XT DNA Sample Preparation Kit (Illumina, San Diego, CA, USA) according to the manufacturer’s protocol [38], unless otherwise stated. Following Illumina’s (Illumina, San Diego, CA, USA) technical note for cluster optimization and the resultant quantity of data each library was normalized for sequencing to 8 pM according to the manufacturer’s protocol. The libraries were quantified using both a Qubit Fluorometer and a 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) and equimolar amounts of DNA were pooled from each sample. Three pools were created with approximately 30 samples per pool, to increase the amount of reads per sample. The three pooled samples were each sequenced at the Hunter College/CTBR Bioinformatics sequencing facility, on an Illumina MiSeq, using the Nano Kit V2 chemistry with 500 cycles and a PhiX spike in of 5%, resulting in 300 bp paired end reads (Illumina, San Diego, CA, USA).

4.4. Bioinformatics Pipeline

Samples were demultiplexed using CASAVA (Illumina, San Diego, CA, USA) allowing no mismatch per barcode. Remaining adapters were trimmed from the reads with Trimmomatic using the provided Truseq adapter sequences [39]. Paired reads were merged using FLASH [40] with a minimum overlap of 250. FastX Toolkit [41] was used to quality filter reads with a minimum of 60% of bases ≥Q30 and a minimum length of 150 bp and was used to trim PCR primer sequences. Trimmed sequences were then clustered at 97% sequence identity using the cd-hit-est algorithm available from CD-HIT [42]. All resulting sequence clusters were screened for chimeras using cd-hit-dup in CD-HIT [42].

4.5. Reads Identification

Obtained sequences were BLASTed against nr database using an e-value cutoff of 1 × 10−5 to discard eventual predator sequences and contaminations. For each sequenced specimen, retained prey sequences were added to a database of ~400 16S rDNA sequences retrieved from NCBI coming from a range of polychaete and sipuncula families. Sequences were aligned using MAFFT [43] and maximum likelihood trees were produced using RAxML [44,45] on the CIPRES portal [46]. Due to the incompleteness of the 16S rDNA polychaetes NCBI database at the species level, we considered it more reliable to identify prey sequences at genus level. All sequences included in each phylogenetic clade were manually inspected and a representative one was chosen and identified with the nearest genus in the tree. Because of the inhomogeneity of our sequencings, we only considered our data to represent the presence or absence of annelid genera. Gut content raw reads (PRJNA662924) and 16S rDNA isolated prey sequences (MW007415–MW007668) are available at GenBank (Supplemental Table S2).

4.6. Diet Data Analysis

For each prey genus detected, we reported the number of predator specimens and species (absolute number and percentage frequency) and the sampling localities, the presence/absence of the venom apparatus, and the phylogenetic clade of the predators.
For estimating the prey width, the Shannon’s diversity index (H’) was calculated with Past3 [47] at the specimen and species level using prey presence/absence data. The Jaccard similarity coefficient, calculated on the same data with Past3 [47], was used to graphically compare predator diets and the results were displayed as heat matrices produced with R [48]. Shannon’s diversity values were then categorized by species, presence/absence of venom apparatus, phylogenetic clade, collecting macro area (Florida, UAE, and Papua New Guinea), and shell length, and obtained distributions were compared using the Kruskal Wallis non-parametric test calculated with Past3 [47]. A randomization test was also conducted to further corroborate the results. Data from all specimens were randomized by column 100 times, the differences in Shannon indices were recalculated, and tested for significance, using R [48]. Using our in-house Terebridae database, we recovered the maximum, minimum, and average collection depth of each species as well as the depth range, and we plotted these data against the species Shannon’s index with Past3 [47].
For investigating correlations between terebrid venom complexity and diet diversity, the transcriptomes from six terebrid specimens (Hastula hectica KVG_110, Hastula matheroniana KVG_193, Terebra subulata KVG_46, Terebra guttata KVG_260, Terebra argus KVG_265, and Terebra straminea KVG_229) were compared with the results of their gut content eDNA. In particular, the number of total putative terebrid toxins and the number of cysteine frameworks obtained following the bioinformatic pipeline from [49] were plotted against the respective average Shannon’s diversity index, using Past3 [47].

Supplementary Materials

The following are available online at https://0-www-mdpi-com.brum.beds.ac.uk/2072-6651/13/2/108/s1. Figure S1: Correlation of terebrid diet and phylogenetic clade and shell length; Table S1: Terebridae vouchers number, foregut type, venom apparatus type, phylogenetic clade, shell length, collection locality, and prey genera found in each specimen; Table S2: Prey sequences number, ID, and GenBank accession number. ID and voucher number of the terebrid from which it was isolated.

Author Contributions

Conceptualization, J.G. and M.H.; Investigation, J.G., G.F. and E.S.L.; Resources, M.H.; Software, J.G.; Visualization, J.G., G.F. and E.S.L.; Writing—original draft, J.G., G.F., E.S.L. and M.H.; Writing—review & editing, J.G., G.F., E.S.L. and M.H. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by Camille & Henry Dreyfus Foundation and the National Institutes of Health (NIH-NIMHD grant 8-G-12-MD007599).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Publicly available datasets were analyzed in this study. This data can be found here: GenBank https://0-www-ncbi-nlm-nih-gov.brum.beds.ac.uk/genbank/ (accessed on 29 January 2021), with accession number: PRJNA662924, MW013507, MW007415–MW007668.

Acknowledgments

M.H. acknowledges research grants by the Camille & Henry Dreyfus Foundation and the National Institutes of Health (NIH-NIMHD grant 8-G-12-MD007599). J.G. acknowledges support from the CUNY Graduate Center Graduate Fellowship program. G.F. acknowledges support from the Regione Lazio Italian fellowship “Torno Subito”. E.S.L. was funded by fellowships from the Undergraduate Research Initiative program and McNulty Scholars program at Hunter College. The authors thank Philippe Bouchet, Virginie Heros, Philippe Maestrati and others involved in obtaining the Papua New Guinea material in this article, which originates from shore-based expeditions and deep-sea cruises (PAPUA NIUGINI), conducted, respectively, by Muséum National d’Histoire Naturelle (MNHN) and Pro-Natura International (PNI) as part of the Our Planet Reviewed program and/or by MNHN, Institut de Recherche pour le Développement (IRD) and other partners as part of the Tropical Deep-Sea Benthos program. In-country scientific partners included the University of Papua New Guinea (UPNG); National Fisheries College (Kavieng, Papua New Guinea). Funders and sponsors included the Total Foundation, Stavros Niarchos Foundation, European Regional Development Fund (ERDF), the French Ministry of Foreign Affairs, and Fonds Pacifique. We acknowledge Maria Vittoria Modica (Stazione Zoologica Anton Dohrn) for advice on the gut content DNA extraction protocol and Nicolas Puillandre (MNHN) for modifying the reverse primer used to amplify prey DNA.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Schendel, V.; Rash, L.D.; Jenner, R.A.; Undheim, E.A. The Diversity of Venom: The Importance of Behavior and Venom System Morphology in Understanding Its Ecology and Evolution. Toxins 2019, 11, 666. [Google Scholar] [CrossRef] [PubMed]
  2. Pekár, S.; Toft, S.; Hrusková, M.; Mayntz, D. Dietary and prey-capture adaptations by which Zodarion germanicum, an ant-eating spider (Araneae: Zodariidae), specialises on the Formicinae. Naturwissenschaften 2008, 95, 233–239. [Google Scholar] [PubMed]
  3. Pucca, M.B.; Amorim, F.G.; Cerni, F.A.; Bordon, K.D.C.F.; Cardoso, I.A.; Anjolette, F.A.P.; Arantes, E.C. Influence of post-starvation extraction time and prey-specific diet in Tityus serrulatus scorpion venom composition and hyaluronidase activity. Toxicon 2014, 90, 326–336. [Google Scholar] [CrossRef] [PubMed]
  4. Daltry, J.C.; Wüster, W.; Thorpe, R.S. Diet and snake venom evolution. Nature 1996, 379, 537–540. [Google Scholar] [CrossRef] [PubMed]
  5. Remigio, E.A.; Duda, T.F., Jr. Evolution of ecological specialization and venom of a predatory marine gastropod. Mol. Ecol. 2008, 17, 1156–1162. [Google Scholar] [CrossRef] [PubMed]
  6. Lyons, K.; Dugon, M.M.; Healy, K. Diet Breadth Mediates the Prey Specificity of Venom Potency in Snakes. Toxins 2020, 12, 74. [Google Scholar] [CrossRef] [PubMed]
  7. Phuong, M.A.; Mahardika, G.N.; Alfaro, M.E. Dietary breadth is positively correlated with venom complexity in cone snails. BMC Genom. 2016, 17, 401. [Google Scholar] [CrossRef]
  8. Vonk, F.J.; Jackson, K.; Doley, R.; Madaras, F.; Mirtschin, P.J.; Vidal, N. Snake venom: From fieldwork to the clinic: Recent insights into snake biology, together with new technology allowing high-throughput screening of venom, bring new hope for drug discovery. Bioessays 2011, 33, 269–279. [Google Scholar] [CrossRef]
  9. Casewell, N.R.; Wüster, W.; Vonk, F.J.; Harrison, R.A.; Fry, B.G. Complex cocktails: The evolutionary novelty of venoms. Trends Ecol. Evol. 2013, 28, 219–229. [Google Scholar] [CrossRef]
  10. Arbuckle, K. From molecules to macroevolution: Venom as a model system for evolutionary biology across levels of life. Toxicon X 2020, 6, 100034. [Google Scholar] [CrossRef]
  11. Li, M.; Fry, B.G.; Kini, R.M. Eggs-only diet: Its implications for the toxin profile changes and ecology of the marbled sea snake (Aipysurus eydouxii). J. Mol. Evol. 2005, 60, 81–89. [Google Scholar] [CrossRef] [PubMed]
  12. Pahari, S.; Bickford, D.; Fry, B.G.; Kini, R.M. Expression pattern of three-finger toxin and phospholipase A2 genes in the venom glands of two sea snakes, Lapemis curtus and Acalyptophis peronii: Comparison of evolution of these toxins in land snakes, sea kraits and sea snakes. BMC Evol. Biol. 2007, 7, 175. [Google Scholar] [CrossRef]
  13. Da Silva, N.J.; Aird, S.D. Prey specificity, comparative lethality and compositional differences of coral snake venoms. Comp. Biochem. Physiol. Part C Toxicol. Pharmacol. 2001, 128, 425–456. [Google Scholar] [CrossRef]
  14. Richards, D.P.; Barlow, A.; Wüster, W. Venom lethality and diet: Differential responses of natural prey and model organisms to the venom of the saw-scaled vipers (Echis). Toxicon 2012, 59, 110–116. [Google Scholar] [CrossRef] [PubMed]
  15. Fry, B.G.; Vidal, N.; Norman, J.A.; Vonk, F.J.; Scheib, H.; Ramjan, S.R.; Kuruppu, S.; Fung, K.; Hedges, S.B.; Richardson, M.K.; et al. Early evolution of the venom system in lizards and snakes. Nature 2006, 439, 584–588. [Google Scholar] [CrossRef] [PubMed]
  16. Fry, B.G.; Casewell, N.R.; Wüster, W.; Vidal, N.; Young, B.; Jackson, T.N. The structural and functional diversification of the Toxicofera reptile venom system. Toxicon 2012, 60, 434–448. [Google Scholar] [CrossRef]
  17. Modica, M.V.; Gorson, J.; Fedosov, A.E.; Malcolm, G.; Terryn, Y.; Puillandre, N.; Holford, M. Macroevolutionary Analyses Suggest That Environmental Factors, Not Venom Apparatus, Play Key Role in Terebridae Marine Snail Diversification. Syst. Biol. 2020, 69, 413–430. [Google Scholar] [CrossRef]
  18. Holford, M.; Puillandre, N.; Terryn, Y.; Cruaud, C.; Olivera, B.; Bouchet, P. Evolution of the Toxoglossa venom apparatus as inferred by molecular phylogeny of the Terebridae. Mol. Biol. Evol. 2009, 26, 15–25. [Google Scholar] [CrossRef]
  19. Castelin, M.; Puillandre, N.; Kantor, Y.I.; Modica, M.V.; Terryn, Y.; Cruaud, C.; Bouchet, P.; Holford, M. Macroevolution of venom apparatus innovations in auger snails (Gastropoda; Conoidea; Terebridae). Mol. Phylogenet. Evol. 2012, 64, 21–44. [Google Scholar] [CrossRef]
  20. Gorson, J.; Holford, M. Small Packages, Big Returns: Uncovering the Venom Diversity of Small Invertebrate Conoidean Snails. Integr. Comp. Biol. 2016, 56, 962–972. [Google Scholar]
  21. Gorson, J.; Ramrattan, G.; Verdes, A.; Wright, E.M.; Kantor, Y.; Rajaram Srinivasan, R. Molecular Diversity and Gene Evolution of the Venom Arsenal of Terebridae Predatory Marine Snails. Genome Biol. Evol. 2015, 7, 1761–1778. [Google Scholar] [CrossRef] [PubMed]
  22. Imperial, J.S.; Kantor, Y.; Watkins, M.; Heralde, F.M., III; Stevenson, B.; Chen, P.; Hansson, K.; Stenflo, J.; Ownby, J.P.; Bouchet, P.; et al. Venomous auger snail Hastula (Impages) hectica (Linnaeus, 1758): Molecular phylogeny, foregut anatomy and comparative toxinology. J. Exp. Zool. B Mol. Dev. Evol. 2007, 308, 744–756. [Google Scholar] [CrossRef] [PubMed]
  23. Anand, P.; Grigoryan, A.; Bhuiyan, M.H.; Ueberheide, B.; Russell, V.; Quinoñez, J.; Moy, P.; Chait, B.T.; Poget, S.F.; Holford, M.; et al. Sample Limited Characterization of a Novel Disulfide-Rich Venom Peptide Toxin from Terebrid Marine Snail Terebra variegata. PLoS ONE 2014, 9, e94122. [Google Scholar] [CrossRef] [PubMed]
  24. Chang, D.; Olenzek, A.M.; Duda, T.F., Jr. Effects of geographical heterogeneity in species interactions on the evolution of venom genes. Proc. Biol. Sci. 2015, 282, 201419874. [Google Scholar] [CrossRef]
  25. Taylor, J.D.; Miller, J.A. A new type of gastropod proboscis: The foregut of Hastula bacillus (Gastropoda: Terebridae). J. Zool. 1990, 220, 603–617. [Google Scholar] [CrossRef]
  26. Kantor, Y.I.; Fedosov, A.E.; Marin, I.N. An Unusually High Abundance and Diversity of the Terebridae (Gastropoda: Conoidea) in Nha Trang Bay, Vietnam. Zool. Stud. 2012, 51, 663–670. [Google Scholar]
  27. Barlow, A.; Pook, C.E.; Harrison, R.A.; Wüster, W. Coevolution of diet and prey-specific venom activity supports the role of selection in snake venom evolution. Proc. Biol. Sci. 2009, 276, 2443–2449. [Google Scholar]
  28. Fry, B.G.; Grigoryan, A.; Bhuiyan, M.H.; Ueberheide, B.; Russell, V.; Quinoñez, J.; Vidal, N.; Poelmann, R.E.; Norman, J.A. Evolution of an arsenal: Structural and functional diversification of the venom system in the advanced snakes (Caenophidia). Mol. Cell. Proteom. 2008, 7, 215–246. [Google Scholar] [CrossRef]
  29. Kordis, D.; Gubensek, F. Adaptive evolution of animal toxin multigene families. Gene 2000, 261, 43–52. [Google Scholar] [CrossRef]
  30. Phuong, M.A.; Mahardika, G.N. Targeted sequencing of venom genes from cone snail genomes reveals coupling between dietary breadth and conotoxin diversity. Mol. Biol. Evol. 2018, 35, 1210–1224. [Google Scholar] [CrossRef]
  31. Van Valen, L. Morphological Variation and Width of Ecological Niche. Am. Nat. 1965, 99, 377–390. [Google Scholar] [CrossRef]
  32. Mans, B.J.; Louw, A.I.; Neitz AW, H. The major tick salivary gland proteins and toxins from the soft tick, Ornithodoros savignyi, are part of the tick Lipocalin family: Implications for the origins of tick toxicoses. Mol. Biol. Evol. 2003, 20, 1158–1167. [Google Scholar] [CrossRef] [PubMed]
  33. Biggs, J.S.; Olivera, B.M.; Kantor, Y.I. Alpha-conopeptides specifically expressed in the salivary gland of Conus pulicarius. Toxicon 2008, 52, 101–105. [Google Scholar] [CrossRef] [PubMed]
  34. Min, G.-S.; Sarkar, I.N.; Siddall, M.E. Salivary transcriptome of the North American medicinal leech, Macrobdella decora. J. Parasitol. 2010, 96, 1211–1221. [Google Scholar] [CrossRef] [PubMed]
  35. Stafford-Banks, C.A.; Rotenberg, D.; Johnson, B.R.; Whitfield, A.E.; Ullman, D.E. Analysis of the Salivary Gland Transcriptome of Frankliniella occidentalis. PLoS ONE 2014, 9, e94447. [Google Scholar] [CrossRef]
  36. Bose, U.; Wang, T.; Zhao, M.; Motti, C.A.; Hall, M.R.; Cummins, S.F. Multiomics analysis of the giant triton snail salivary gland, a crown-of-thorns starfish predator. Sci. Rep. 2017, 7, 6000. [Google Scholar] [CrossRef]
  37. Duda, T.F., Jr.; Chang, D.; Lewis, B.D.; Lee, T. Geographic Variation in Venom Allelic Composition and Diets of the Widespread Predatory Marine Gastropod Conus ebraeus. PLoS ONE 2009, 4, e6245. [Google Scholar] [CrossRef]
  38. Illumina Nextera (R) XT DNA Library Preparation Guide. 2014.
  39. Bolger, A.M.; Lohse, M.; Usadel, B. Trimmomatic: A flexible trimmer for Illumina sequence data. Bioinformatics 2014, 30, 2114–2120. [Google Scholar] [CrossRef]
  40. Magoč, T.; Salzberg, S.L. FLASH: Fast length adjustment of short reads to improve genome assemblies. Bioinformatics 2011, 27, 2957–2963. [Google Scholar] [CrossRef]
  41. Lab, H. FASTX Toolkit; Cold Spring Harbor Laboratory: New York, NY, USA, 2014. [Google Scholar]
  42. Li, W.; Godzik, A. Cd-hit: A fast program for clustering and comparing large sets of protein or nucleotide sequences. Bioinformatics 2006, 22, 1658–1659. [Google Scholar] [CrossRef]
  43. Katoh, K.; Misawa, K.; Kuma, K.-I.; Miyata, T. MAFFT: A novel method for rapid multiple sequence alignment based on fast Fourier transform. Nucleic Acids Res. 2002, 30, 3059–3066. [Google Scholar] [CrossRef] [PubMed]
  44. Stamatakis, A. RAxML-VI-HPC: Maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models. Bioinformatics 2006, 22, 2688–2690. [Google Scholar] [CrossRef] [PubMed]
  45. Stamatakis, A. RAxML version 8: A tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics 2014, 30, 1312–1313. [Google Scholar] [CrossRef]
  46. Miller, M.A.; Pfeiffer, W.; Schwartz, T. Creating the CIPRES Science Gateway for inference of large phylogenetic trees. Gatew. Comput. Environ. Workshop 2010. [Google Scholar] [CrossRef]
  47. Hammer, Ø.; Harper, D.A.T.; Ryan, P.D. PAST: Paleontological statistics software package for education and data analysis. Palaeontol. Electron. 2001, 4, 9. [Google Scholar]
  48. R Core Team. R: A Language and Environment for Statistical Computing; R Foundation for Statistical Computing: Vienna, Austria, 2017; Available online: https://www.r-project.org/ (accessed on 29 January 2021).
  49. Eriksson, A.; Anand, P.; Gorson, J.; Grijuc, C.; Hadelia, E.; Stewart, J.C.; Holford, M.; Claridge-Chang, A. Using Drosophila behavioral assays to characterize terebrid venom-peptide bioactivity. Sci. Rep. 2018, 8, 15276. [Google Scholar] [CrossRef]
Figure 1. Representative venomous terebrid snails and their annelid prey identified in this study. (1) Hastula hectica. (2) Myurella nebulosa. (3) Terebra subulata. (4) Terebra argus. (5) Myruellopsis undulata. (6) Terebra guttata. (7) Scolelepis. (8) Neanthes. (9) Pygospio. (10) Odontosyllis. (11) Syllis. (12) Worm pulled from the mouth of Terebra guttata (6), identified as Scolelepis. Photo credits available online: https://www.aphotomarine.com (accessed on 29 January 2021).
Figure 1. Representative venomous terebrid snails and their annelid prey identified in this study. (1) Hastula hectica. (2) Myurella nebulosa. (3) Terebra subulata. (4) Terebra argus. (5) Myruellopsis undulata. (6) Terebra guttata. (7) Scolelepis. (8) Neanthes. (9) Pygospio. (10) Odontosyllis. (11) Syllis. (12) Worm pulled from the mouth of Terebra guttata (6), identified as Scolelepis. Photo credits available online: https://www.aphotomarine.com (accessed on 29 January 2021).
Toxins 13 00108 g001
Figure 2. Eighteen different annelid prey genera found in gut of terebrid snails. Shown is the phylogeny of terebrid species used in this study and the annelid worm meal identified in the gut of terebrid specimens. Annelid meals were identified for terebrid species with a venom gland and for those without a venom gland. Terebrid species with a traditional venom gland are highlighted in red in the tree. Smaller pie charts represent the diet of each adjacent taxon. Larger pie charts represent the aggregate diet for terebrids in the same clade. The phylogeny is adapted from the work of Modica et al., 2020, the most up-to-date and first dated phylogeny of the Terebridae family [17].
Figure 2. Eighteen different annelid prey genera found in gut of terebrid snails. Shown is the phylogeny of terebrid species used in this study and the annelid worm meal identified in the gut of terebrid specimens. Annelid meals were identified for terebrid species with a venom gland and for those without a venom gland. Terebrid species with a traditional venom gland are highlighted in red in the tree. Smaller pie charts represent the diet of each adjacent taxon. Larger pie charts represent the aggregate diet for terebrids in the same clade. The phylogeny is adapted from the work of Modica et al., 2020, the most up-to-date and first dated phylogeny of the Terebridae family [17].
Toxins 13 00108 g002
Figure 3. Terebrid snails without a venom gland (VG) have a larger diversity in prey diet. (A) Diet diversity in terebrids with and without a VG. (B) Histogram showing the frequency of Shannon index differences obtained from Kruskal Wallis non-parametric test on data categorized by venom apparatus (randomized data in blue and study data in red).
Figure 3. Terebrid snails without a venom gland (VG) have a larger diversity in prey diet. (A) Diet diversity in terebrids with and without a VG. (B) Histogram showing the frequency of Shannon index differences obtained from Kruskal Wallis non-parametric test on data categorized by venom apparatus (randomized data in blue and study data in red).
Toxins 13 00108 g003
Figure 4. No direct correlation between terebrid diet and venom complexity. Diet breadth, calculated as the Shannon’s diversity index (H’), was plotted against the number of putative toxins (A) and with the number of cysteine frameworks (B) in six terebrids. No correlation was detected. Analyses were performed with Past3.
Figure 4. No direct correlation between terebrid diet and venom complexity. Diet breadth, calculated as the Shannon’s diversity index (H’), was plotted against the number of putative toxins (A) and with the number of cysteine frameworks (B) in six terebrids. No correlation was detected. Analyses were performed with Past3.
Toxins 13 00108 g004
Figure 5. Terebrid species with different diet produce similar venom toxins. Heat map showing the interspecific variation of relative cysteine frameworks in terebrid toxins. Adjacent pie charts represent the identified diet for each specimen. Most of the specimens show similar relative percentages of each cysteine frameworks but H. hectica has higher values for frameworks I and XI, and no framework XXII putative teretoxins found in all other terebrid samples use in this study.
Figure 5. Terebrid species with different diet produce similar venom toxins. Heat map showing the interspecific variation of relative cysteine frameworks in terebrid toxins. Adjacent pie charts represent the identified diet for each specimen. Most of the specimens show similar relative percentages of each cysteine frameworks but H. hectica has higher values for frameworks I and XI, and no framework XXII putative teretoxins found in all other terebrid samples use in this study.
Toxins 13 00108 g005
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Back to TopTop