# Redistribution of Ancestral Functions Underlies the Evolution of Venom Production in Marine Predatory Snails

Giulia Zancolli, Maria Vittoria Modica, Nicolas Puillandre, Yuri Kantor, Agneesh Barua, Giulia Campli, Marc Robinson-Rechavi

► To cite this version:

Giulia Zancolli, Maria Vittoria Modica, Nicolas Puillandre, Yuri Kantor, Agneesh Barua, et al.. Redistribution of Ancestral Functions Underlies the Evolution of Venom Production in Marine Predatory Snails. Molecular Biology and Evolution, 2025, 42 (5), pp.msaf095. 10.1093/molbev/msaf095. hal-05214272

HAL Id: hal-05214272

https://hal.science/hal-05214272v1

Submitted on 19 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Redistribution of Ancestral Functions Underlies the Evolution of Venom Production in Marine Predatory Snails

Giulia Zancolli ID,<fnref n="1" />,<fnref n="2" />,* Maria Vittoria Modica ID,<fnref n="3" /> Nicolas Puillandre ID,<fnref n="4" /> Yuri Kantor ID,<fnref n="5" />

Agneesh Barua ID,<fnref n="1" />,<fnref n="2" /> Giulia Campli ID,<fnref n="1" />,<fnref n="2" /> Marc Robinson-Rechavi ID,<fnref n="1" />,<fnref n="2" />

## Abstract

Venom-secreting glands are highly specialized organs evolved throughout the animal kingdom to synthesize and secrete toxins for predation and defense. Venom is extensively studied for its toxin components and application potential; yet, how animals become venomous remains poorly understood. Venom systems therefore offer a unique opportunity to understand the molecular mechanisms underlying functional innovation. Here, we conducted a multispecies multi-tissue comparative transcriptomics analysis of 12 marine predatory gastropod species, including species with venom glands and species with homologous non-venom-producing glands, to examine how specialized functions evolve through gene expression changes. We found that while the venom gland specialized for the mass production of toxins, its homologous glands retained the ancestral digestive functions. The functional divergence and specialization of the venom gland were achieved through a redistribution of its ancestral digestive functions to other organs, specifically the esophagus. This entailed concerted expression changes and accelerated transcriptome evolution across the entire digestive system. The increase in venom gland secretory capacity was achieved through the modulation of an ancient secretory machinery, particularly genes involved in endoplasmic reticulum stress and unfolded protein response. This study shifts the focus from the well-explored evolution of toxins to the lesser-known evolution of the organ and mechanisms responsible for venom production. As such, it contributes to elucidating the molecular mechanisms underlying organ evolution at a fine evolutionary scale, highlighting the specific events that lead to functional divergence.

**Keywords:** Neogastropoda, transcriptome evolution, ancestral reconstruction, toxin, venom, UPR, ER stress, secretion

## Introduction

Across diverse branches of the animal kingdom, organisms have independently evolved the ability to produce and deliver venom, a cocktail of bioactive toxin molecules. In many venomous animals, these toxins are synthesized in specialized exocrine organs known as venom glands (Casewell et al. 2013; Schendel et al. 2019). While significant research has focused on the molecular evolution of toxins and venom composition, the molecular mechanisms underlying the evolution of venom-producing organs and their specialized function remains poorly understood (Zancolli and Casewell 2020). Recent genomics studies in snakes, for instance, have highlighted how regulatory networks are co-opted to drive high toxin gene expression in venom glands, notably through the involvement of trans-regulatory factors from the extracellular signal-regulated kinase (ERK) and unfolded protein response (UPR) pathways (Perry et al. 2020, 2022; Westfall et al. 2023). Moreover, the upregulation of UPR and endoplasmic reticulum (ER) stress pathways in the venom glands of several distinct venomous taxa suggests that similar molecular solutions may have convergently evolved across venomous lineages to support the high secretory demands of toxin

production (Perry et al. 2020; Barua and Mikheyev 2021; Zancolli et al. 2022). However, a broader understanding of how venom glands become highly specialized and optimized for the efficient mass production of toxins is still lacking, particularly outside of snake models.

Investigating homologous structures with divergent functions provides valuable insights into the genetic basis of functional innovation, as demonstrated in large-scale evolutionary studies on vertebrate limbs (Zuniga 2015) and feathers (Benton et al. 2019). Among non-vertebrates, the mid-esophageal glands of marine predatory snails in the subclass Caenogastropoda (Fig. 1) offer an excellent system to explore, at a finer evolutionary scale, the specific events that led to the emergence of new physiological functions.

The digestive system of carnivorous snails is highly diverse and characterized by the presence of glandular structures in the middle region of the esophagus (Lobo-da-Cunha 2019). The development and homology of these structures, referred to as mid-esophageal glands, have been extensively studied by Page and others, who demonstrated that in the larvae of various caenogastropods, the hypertrophy of epithelial cells in a ventral zone of the mid-foregut region gives rise to the mid-esophageal glands (Page 2000, 2002, 2012; Page and

<footnote n="1">¹Department of Ecology and Evolution, University of Lausanne, 1015 Lausanne, Switzerland</footnote>

<footnote n="2">²Evolutionary Bioinformatics, Swiss Institute of Bioinformatics, Lausanne, Switzerland</footnote>

<footnote n="3">³Department of Biology and Evolution of Marine Organisms, Stazione Zoologica Anton Dohrn, 00198 Roma, Italy</footnote>

<footnote n="4">⁴Institut Systématique Evolution Biodiversité (ISYEB), Muséum National d'Histoire Naturelle, CNRS, Sorbonne Université, EPHE, Université des Antilles, 75005 Paris, France</footnote>

<footnote n="5">⁵Severtsov Institute of Ecology and Evolution, Russian Academy of Sciences, 119034 Moscow, Russian Federation</footnote>

<footnote>*Corresponding author: E-mail: giulia.zancolli@gmail.com.</footnote>

<footnote>Associate editor: Ilya Ruvinsky</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Species and tissues investigated. Overview of the anatomy of a marine predatory gastropod with the sampled tissues and the phylogeny of the species used in this study based on Fedosov et al. (2024). The abbreviations used for the tissues and species used in other figures are shown. The anatomy drawing is modified after (Modica and Holford 2010), while the schematic of the foregut apparatus of caenogastropods with their mid-esophageal glands was modified after (Page 2011).</figcaption></figure>

Hookham 2017). In some caenogastropods, like those in the family Naticidae, the mid-esophageal gland consists of a di- lated section of the esophageal wall involved in secreting digestive enzymes and mucus (Reid and Friesen 1980; Lobo-da-Cunha 2019). By contrast, in many caenogastropods of the order Neogastropoda, this glandular section of the esophagus has evolved into a distinct organ known as the gland of Leiblein, which connects to the esophagus via a duct (Modica and Holford 2010; Lobo-da-Cunha 2019). Ultrastructural observations in species from the families Muricidae and a Nassariidae species suggest that the gland of Leiblein plays a role in food processing, particularly nutri- ent absorption and storage (Andrews and Thorogood 2005). In the superfamily Conoidea, which comprises 18 families (Abdelkrim et al. 2018) including the well-known cone snails (family Conidae) (Nguyen et al. 2023), the esophageal gland has undergone further modifications into a long, convoluted duct that secretes a mixture of hundreds of primarily neurotoxic peptides, mainly known as “conotoxins” (Modica and Holford 2010; Abdelkrim et al. 2018; Morales Duque et al. 2019). The venom gland is attached to a large muscular bulb that contracts to push venom through the duct and into the buccal cavity, where a modified radula injects venom into prey or predators (Modica and Holford 2010; Safavi-Hemami et al. 2010; Salisbury et al. 2010). Interestingly, some neogastropods, such as those in the family Mitridae (Harasewych 2009) and Terebridae (Holford et al. 2009), have entirely lost the mid-esophageal gland.

The evolution of venom production from digestive-related
functions in these snails represents a remarkable example of
functional specialization, shifting the gland's secretion targets
from internal digestive roles (e.g. lysosomal activity) (Andrews
and Thorogood 2005) to external roles that affect other or-
ganisms (e.g. receptor binding in prey) (Morales Duque et al.
2019). This transition likely provided significant adaptive ad-
vantages, enabling cone snails to diversify their diet to include
fast-moving organisms, such as the recent adaptation to fish
hunting, while also offering a defense against powerful preda-
tors (Olivera et al. 2012). However, the processes through
which an organ originally dedicated to digestive functions

transformed into a specialized toxin-producing factory remain
unclear.

Here, we analyze gene expression data from the mid-esophageal glands and other tissues of 12 marine caenogastropod species to investigate the link between transcriptome evolution and functional divergence. By studying the molecular underpinnings of venom gland evolution, we aim to uncover how specific genes, pathways, and regulatory networks drive organ specialization to meet the distinct physiological demands of venom production. Our study addresses the following key questions: (i) Do mid-esophageal glands across species share similar gene expression profiles, given their common origin? (ii) Considering the venom gland's unique function, does its transcriptome evolve more rapidly compared to the other mid-esophageal glands? (iii) Which gene expression changes led to the evolution of toxin production in the mid-esophageal gland?

To answer these questions, we first conducted species-level
analyses to characterize sets of overexpressed genes and to de-
lineate the functional specialization of each gland type. We
then performed between-species comparisons to explore the
gene expression dynamics that led to the evolution of venom
production. Our findings reveal that the venom gland has a
markedly distinct gene expression profile compared to its
homologous organs, which are more similar to each other.
Genes encoding secreted proteins in venom glands are ex-
pressed at exceptionally high levels, far exceeding even those
of other secretory organs. This specialization for toxin secretion
was achieved through modulation of a conserved secretory ma-
chinery, while ancestral digestive functions were redistributed
to other organs. This shift involved high evolutionary rates
not only in the venom gland itself but across the entire digestive
system, suggesting concerted changes that underscore the adap-
tive flexibility of these organisms.

Results

Summary of Species and Tissues Investigated

We sampled ten species of Neogastropoda and two outgroup
species within the same subclass Caenogastropoda (Fig. 1).
The two outgroup species, the Cymatidae *Monoplex pilearis* (Linnaeus, 1758) and the Naticidae *Polinices mammilla* (Linnaeus, 1758), possess a simple esophageal gland (OEG) attached to the esophagus, representing the ancestral state. Among the Neogastropoda, two species from the Muricidae family, *Nassa serta* (Bruguière, 1789) and *Murex tenuirostrum* (Lamarck, 1822), have a gland of Leiblein (LEG), with *M. tenuirostrum* also possessing dorsal glandular folds on the mid-esophagus, often referred to as glande framboisée (Lobo-da-Cunha 2019). Additionally, we sampled the gland of Leiblein in the Olividae *Oliva amethystina* (Röding, 1798), the Vasidae *Vasum turbinellus* (Linnaeus, 1758), and the Nassaridae *Phos senticosus* (Linnaeus, 1758). The venomous species, possessing a venom gland (VG), belonged to the family Conidae, including *Conus imperialis* (Linnaeus, 1758), *Conus marmoreus* (Linnaeus, 1758), *Conus virgo* (Linnaeus, 1758), and *Conus quercinus* (Lightfoot, 1786). We also collected tissue samples from the Mitridae *Mitra mitra* (Linnaeus, 1758), which lacks a mid-esophageal gland (Harasewych 2009). Besides the glands, we sampled either the foot or the columellar muscle, esophagus, salivary glands, dorsal glandular folds in *M. tenuirostrum*, and muscular venom bulb in cone snails.

## Sequencing and De Novo Assembly Statistics

A total of 150 libraries were sequenced, yielding an average of 45 million reads per library across 12 species, with an average of three biological replicates per species (supplementary dataset S1, Supplementary Material online). We generated de novo assemblies by pooling all libraries within each species and processing them through our quality-filtering pipeline (see Materials and Methods). After filtering, the assemblies contained between 26,215 and 59,431 annotated, nonredundant transcripts, with an average of 40,192 sequences per assembly (supplementary dataset S2, Supplementary Material online). The completeness of these assemblies ranged from 82% to 93%, with an average of 89%.

Prior to data analysis, we evaluated the quality of read alignment. Samples with low mapping rates were excluded. Furthermore, we confirmed that samples of the same tissue type clustered together in a Principal Component Analysis (PCA) plot (supplementary fig. S1, Supplementary Material online). After this quality filtering step, we retained a total of 140 samples.

## Within-Species Analysis

### Characterization of Tissue-Specific Gene Sets

To better understand the functional specialization of the mid-esophageal glands, we identified sets of genes that were over-expressed in each organ. Genes were classified as tissue-specific if their expression in a tissue was at least twice that of the second most highly expressed tissue. On average, 14% (range: 10% to 19%) of genes were tissue-specific, with 43% (range: 38% to 51%) of these genes exclusively expressed in one tissue. Tissue specificity was consistent across organs and species, although we observed a higher number of tissue-specific genes in the venom glands compared to the glands of Leiblein and esophageal glands (Fig. 2a). The smallest set observed was in the esophagus of *M. tenuirostrum* (N=318). Interestingly, this species possesses dorsal glandular folds on the esophagus (Fig. 1). The tissue specificity of the gene sets was further validated by differential expression

analysis using a likelihood ratio test across all tissues (supplementary text S1, Supplementary Material online).

Gene Ontology (GO) enrichment analysis of the tissue-specific gene sets revealed marked differences between the mid-esophageal glands (Fig. 2b), as well as lineage-specific patterns (supplementary text S2 and figs. S2 to S10, Supplementary Material online). In OEG-specific gene sets, we observed enrichment for terms related to intracellular trafficking and communication (e.g. “AP-5 adaptor complex”) and transmembrane transport. Additionally, terms related to detoxification, homeostasis, and response to external stimuli (e.g. “response to inorganic substance” and “superoxide dismutase activity”) were enriched (supplementary figs. S2 to S4, Supplementary Material online). In contrast, gene sets specific to the gland of Leiblein across multiple species were enriched for terms related to intracellular digestion, particularly protein digestion and absorption, with several lysosome-related terms (e.g. “endopeptidase activity” and “WASH complex”), aligning with previous ultrastructural studies (Andrews and Thorogood 2005). Terms related to iron homeostasis and metabolism were also enriched. However, *O. amethystina* showed some unique enrichments, including extracellular rather than intracellular compartment terms, reflecting possible lineage-specific adaptations linked to diet (supplementary figs. S5 to S7, Supplementary Material online). In venom glands, GO terms associated with toxin and neurotoxic activity (e.g. “host cell postsynaptic membrane”) were over-represented in all species, as expected (Fig. 2b). Additionally, we found terms related to protein synthesis (e.g. “signal peptide processing” and “Golgi organisation”) and to post-translational modifications (e.g. “hydroxylation” and “dissulfide oxidoreductase activity”). Notably, terms related to ER stress and UPR were also enriched (e.g. “response to unfolded proteins” and “response to ER stress”) (supplementary figs. S8 to S10, Supplementary Material online).

In summary, while the nonvenomous mid-esophageal glands are primarily involved in digestion, homeostasis, absorption, and storage, the homologous venom gland has specialized into a factory for toxin synthesis and secretion. Despite these differences, all three gland types shared enriched terms related to iron homeostasis and metabolism (e.g. “iron binding”) and methylation (“methionine adenosyltransferase activity”, “betaine-homocysteine S-methyltransferase activity”). Additionally, we found terms related to hormone response across all three gland types, including “thyroid hormone generation” in OEG- and LEG-specific gene sets, and “response to thyroglobulin triiodothyronine” in VG-specific sets.

### Characterization of the Glands' Secretomes

Given the fundamental role of secretion in the evolution of the venom gland, we analyzed the “secretome” by assessing the diversity and expression levels of genes predicted to have a signal peptide with SignalP (Teufel et al. 2022).

Gene expression of secretory proteins was relatively consistent across tissue types, although the number of expressed genes was higher in the venom gland (mean N=565) compared to the gland of Leiblein (t-test: mean N=354, t=2, df=5.8, P=0.04) and the esophageal gland (mean N=379, t=1, df=1,3, P=0.23), although the latter was not statistically significant (supplementary fig. S11a, Supplementary Material online). When focusing on tissue-specific genes, the venom gland secretome was also significantly more diverse (mean N=183) than the other glands (gland of Leiblein: mean=N=55; t=2.2, df=3.3, P=0.05; esophageal gland:
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Overview of tissue-specific gene sets and gland secretomes. a) Number of genes within each tissue-specific set across all organs (x axis) and species (y axis). Species abbreviations as in Fig. 1. The venom glands have a significant higher number of tissue-specific genes (mean N = 1,424) compared to the gland of Leiblein (one-tailed t-test: mean N = 771; t = 3.6, df = 6, P = 0.005) and esophageal gland (mean N = 699; t = 5.6, df = 3, P = 0.005). Missing tissue samples are marked with a “-”. b) Gene Ontology (GO) enrichment results of the OEG-, LEG-, and VG-specific gene sets. c) Mean expression levels as Transcript per Million (TPM) of genes possessing a signal peptide, therefore comprising the secretomes, expressed in the organs (x axis) of each species (y axis). Species abbreviations as in Fig. 1. Missing tissue samples are marked with a “-”. d) GO enrichment results of the esophageal gland, gland of Leiblein, and venom gland secretomes.</figcaption></figure>

mean N = 38; t = 2.5, df = 3, P = 0.04) (supplementary fig. S11b, Supplementary Material online). Interestingly, the salivary glands did not show particularly high diversity despite also being exocrine glands like the venom gland (supplementary fig. S11, Supplementary Material online).

In terms of expression levels (mean TPM) of the secretome, the differences between tissues were more pronounced (Fig. 2c). Gene expression of secreted proteins was generally higher in venom glands, although it varied among cone snail species (Fig. 2c). This variation could reflect genuine lineage-specific differences, technical factors, or differences in the venom replenishment cycle at the time of sampling. However, the latter is unlikely, as individuals were sampled randomly and thus would not all be at the same point in their replenishment cycle. Additionally, all specimens were kept for 1 to 3 d in captivity prior dissection to minimize environmental variation. We also did not observe ingested prey in any stomachs, suggesting that the last feeding—and thus venom expulsion—did not occur close to the time of dissection for

any individual. In nonvenomous species, the highest secretome expression levels were observed in the salivary glands. Interestingly, in the esophagus, lower expression levels were observed in nonvenomous species (mean = 21 TPM) compared to venomous species (one-tailed t-test, mean = 260; t = -3, df = 3, P = 0.02), suggesting higher secretory activity in the latter.

As anticipated, the secretomes of the esophageal gland and gland of Leiblein were enriched in hydrolases and peptidases, enzymes essential for digestive processes (Fig. 2d). Additionally, *O. amethystina* showed enrichment in “toxin activity” (see next section), while *M. pilearis* was enriched in terms related to communication and transport (e.g. “metal ion transport”). Although the salivary glands also expressed hydrolases and peptidases, these enzymes were extracellular, whereas those in the gland of Leiblein were primarily intracellular, consistent with enrichment in cellular compartments like “lysosome” and “organelle lumen”. In contrast, the venom gland secretomes were dominated by toxins released in
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Transcriptome similarity and shared tissue specificity. a) Heatmap of Pearson correlation coefficients between tissues and species. The expression tree was made using neighbor-joining based on the correlation matrix. b) Average number of OEG-, LEG-, and VG-specific OGs shared with other tissue-specific OG sets, where species have been grouped by their gland type (OEG-, LEG-, VG- species and glandless).</figcaption></figure>

the extracellular space and of genes involved in the ER function (Fig. 2d).

## Identification and Characterization of Conotoxins

Cone snail venom is composed primarily of small, disulfide-rich peptides known as conotoxins. The number of distinct putative conotoxin transcripts predicted in the venomous species' assemblies was consistent with previous findings for cone snail venom gland de novo transcriptomes (Abalde et al. 2018; Dutt et al. 2019), with 150 to 250 toxins predicted per species. Notably, conotoxin-like sequences were also identified in other species, with numbers ranging from 33 in *M. mitra* to 88 in *O. amethystina*. However, when restricting the analysis to sequences with both a signal peptide and a predicted conotoxin domain, the numbers were greatly reduced, ranging from 16 to 108 in venomous species, and 7 to 28 in nonvenomous species.

As expected, conotoxins were predominantly expressed in the venom gland (supplementary fig. S12, Supplementary Material online). A few were also highly and specifically expressed in the salivary glands, consistent with previous reports in other cone snail species (Fedosov et al. 2023). In nonvenomous species, predicted conotoxins were expressed at much lower levels and across multiple tissues, although a trend toward higher expression in salivary glands was observed (supplementary figs. S13 and S14, Supplementary Material online).

It is important to note that the prediction tool used, ConoPrec (Kaas et al. 2012), is designed to predict only conoidean toxins. Consequently, toxins found outside Conoidea, such as echotoxin (Shiomi et al. 2002), were not predicted, leading to an underestimation of the toxic potential of gastropod salivary glands (Ponte and Modica 2017). Our focus on conotoxins stems from their critical role as the primary weapon of cone snails, whose massive gene expansion and diversification likely drove the evolution of the venom gland.

## Between-Species Analysis

### Orthogroup Assignment

We assigned 330,491 genes (68% of the total) to 41,720 orthogroups (OGs), with 2,588 OGs shared across all species. For comparative transcriptomics, we created a multispecies expression matrix using these 2,588 OGs. Since many OGs contained multiple genes per species, we selected a representative gene for each OG using two methods: (i) randomly selecting a single transcript's TPM value, and (ii) calculating the mean TPM across all the transcripts within an OG. Both approaches produced similar outcomes, therefore the results presented here are based on the first method. Results from the second approach are provided in the Supplementary Material online.

### Transcriptome Similarity and Shared Tissue Specificity Between Organs and Species

To determine whether homologous glands share similar global gene expression profiles or if their functional specializations align them more closely with nonhomologous organs, we analyzed whole transcriptome similarity patterns by means of correlation matrix and PCA. Overall, samples primarily clustered by tissue type, with the exception of the esophagus and salivary glands samples, which showed more overlap (Fig. 3a, supplementary figs. S15 and S16, Supplementary Material online). The esophageal glands and glands of Leiblein grouped together, while the venom glands clustered with the salivary glands rather than with their homologous counterparts.

Interesting patterns emerged when comparing tissue specificity across species (Fig. 3b). Specifically, we calculated the average from pairwise comparisons of gland-specific OGs between all species. We found that the esophageal glands on average shared more tissue-specific OGs with the glands of Leiblein (mean N = 61, SD = 17.7) and the esophagus of the
<table><caption>**Table 1** Evolutionary rates for the tested models</caption><thead><tr><th>Tissue</th><th>Model</th><th>N. of evolutionary rates (σ<sup>2</sup>)</th><th>Likelihood (-ln L)</th><th>OEG-species σ<sup>2</sup></th><th>LEG-species σ<sup>2</sup></th><th>Glandless-species σ<sup>2</sup></th><th>VG-species σ<sup>2</sup></th></tr></thead><tbody><tr><td rowspan="3">Gland</td><td rowspan="2">Random</td><td>3</td><td>17,109.7</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>1</td><td>17,649.4</td><td>1.02</td><td>1.02</td><td>...</td><td>1.02</td></tr><tr><td>2</td><td>2</td><td>16,525.1</td><td>0.59</td><td>0.59</td><td>...</td><td>1.79</td></tr><tr><td rowspan="3">Salivary glands</td><td rowspan="2">Random</td><td>3</td><td>16,018.7</td><td>0.23</td><td>0.72</td><td>...</td><td>1.80</td></tr><tr><td>3</td><td>17,480</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>1</td><td>1</td><td>17,953.2</td><td>0.93</td><td>0.93</td><td>0.93</td><td>0.93</td></tr><tr><td rowspan="3">Esophagus</td><td rowspan="2">Random</td><td>2</td><td>15,434.8</td><td>0.51</td><td>0.51</td><td>0.51</td><td>2.93</td></tr><tr><td>3</td><td>15,314.9</td><td>0.35</td><td>0.58</td><td>0.35</td><td>2.86</td></tr><tr><td>1</td><td>1</td><td>22,018.2</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td rowspan="3"></td><td rowspan="2">Random</td><td>2</td><td>22,728.1</td><td>1.15</td><td>1.15</td><td>1.15</td><td>1.15</td></tr><tr><td>3</td><td>21,316.9</td><td>0.70</td><td>0.70</td><td>0.70</td><td>2.11</td></tr><tr><td>3</td><td>21,089.5</td><td>0.39</td><td>0.81</td><td>0.39</td><td>2.10</td></tr></tbody></table>

For each organ is reported: the number of evolutionary rates estimated (σ²), the likelihood, and the σ² estimates for each group. The best model (i.e. lowest likelihood value) is highlighted in bold.

glandless *M. mitra* (mean N = 56, SD = 10.6) rather than with each other (mean N = 41, SD = 0). This suggests substantial variation between the two OEG-species, as corroborated by the GO enrichment analysis (supplementary figs. S2 to S4, Supplementary Material online). The glands of Leiblein, in contrast, shared more tissue-specific OGs among themselves (mean N = 103, SD = 45.9), and with the esophagus of *M. mitra* (mean N = 100, SD = 38.2). In contrast, venom glands shared more tissue-specific OGs among themselves (mean N = 111, SD = 20.4) and with the salivary glands of *M. mitra* (mean N = 71, 11.5), while very few with the esophageal glands (mean N = 11, SD = 3.6) and glands of Leiblein (mean N = 16, 7.7). In summary, while the esophageal gland and gland of Leiblein share more similar transcriptomes, the venom gland markedly diverged from its homologous organs and converged toward the other exocrine organ, the salivary glands (Fig. 3b).

For subsequent analyses, we defined OGs as tissue-specific if they were upregulated in the same tissue of at least two species, resulting in 41 OEG-specific OGs, 343 LEG-specific OGs, and 405 VG-specific OGs.

## Rates of Gene Expression Evolution

Given the marked functional divergence of the venom gland, we hypothesized that its transcriptome evolves faster than that of its homologous glands. We tested this hypothesis using CAGEE (Bertram et al. 2023), which employs a bounded Brownian motion model to estimate the most likely value of the evolutionary rate parameter (σ²) consistent with an ultra-metric species tree and with the observed expression values at the tip of the tree. The evolutionary rate is then used to infer genome-wide changes in gene expression along the phylogenetic tree, and to calculate ancestral states for each gene at each node (Bertram et al. 2023).

We ran CAGEE for the mid-esophageal glands, salivary glands, and esophagus, as these tissues were sampled across most species representing all three gland types. We evaluated four different evolutionary models (supplementary fig. S17, Supplementary Material online). The first model estimated a single rate σ². The second model estimated two distinct rates, one for the venomous clade and one for all other species. In the third model, species were grouped by gland type and the rates were estimated for each group separately. The final model also

calculated three rates but assigned them randomly across the phylogeny. Overall, the third model had the best fit (Table 1, supplementary table S1, Supplementary Material online). Notably, the model with the poorest fit was the first one, indicating that some degree of variation in evolutionary rates, even if random, is more consistent with the data than assuming a uniform rate across all lineages. When comparing across lineages, the venom gland showed the highest rate (σ² value), supporting our hypothesis of accelerated evolution in the venom gland relative to the other homologous glands. A similar trend was observed in other organs, with higher evolutionary rates in venomous species than nonvenomous ones. However, when comparing across organs, none of the glands had the highest σ² value, suggesting that other organs also underwent accelerated evolution, even more so than the mid-esophageal glands, supporting the hypothesis of concerted evolution across the entire digestive system.

## Expression Changes of Tissue-Specific Orthogroups Along the Phylogeny

Given the better fit of the third evolutionary model, we used this one for ancestral state reconstruction to assess the number and direction of expression changes at each node of the phylogenetic tree (Fig. 4b). We observed substantial changes at nodes 17 and 18, which lead to the glandless *M. mitra* and the venomous clade (Fig. 4a, supplementary fig. S18, Supplementary Material online). At these nodes, many OEG- and LEG-specific OGs showed a marked decrease in expression in the ancestral gland, coupled with an increase in expression in the esophagus and salivary glands. Conversely, VG-specific OGs underwent upregulation in the gland as well as across other tissues (Fig. 4a).

Among the OEG- and LEG-specific OGs, the strongest downregulation in the gland was a galectin, which simultaneously showed the highest upregulation in the esophagus (Fig. 4d), alongside with a member of the ependymin family. The top three VG-specific OGs that were upregulated in the gland included an integral membrane protein of the DAD family, a disulfide-isomerase, and a selenoprotein F (Fig. 4d), all involved in the ER.

We then utilized the ancestral state reconstruction from CAGEE to calculate tissue specificity at each node, akin to our approach for extant species. The goal was to understand
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Gene expression dynamics across the phylogeny. a) Number of OEG-, LEG, and VG-specific OGs decreasing and increasing their expression levels in the ancestral salivary glands, esophagus, or mid-esophageal gland at each internal node of the gastropod phylogeny. The expression changes were calculated based on gene expression reconstruction at each node of the phylogeny. b) Species phylogeny with the number of the internal nodes. c) Number of OEG-, LEG-, and VG-specific OGs which were tissue-specific in the ancestral mid-esophageal gland, salivary glands, and esophagus at each node of the phylogeny. d) Ancestral reconstruction of gene expression of a galectin (OG2631) and selenoprotein F (OG465) in the mid-esophageal gland and esophagus. The missing tissues are marked with a *.</figcaption></figure>

whether gland-specific OGs were tissue-specific in the ances-
tral lineages. Our findings reveal that some OEG- and
LEG-specific OGs were gland-specific in ancestral lineages un-
til node 18, where they shifted mainly to the esophagus
(Fig. 4c, supplementary fig. S18, Supplementary Material on-
line). The galectin and ependymin mentioned earlier were
among the genes that transitioned to esophagus specificity,
alongside other genes involved in gluconeogenesis. In contrast,
VG-specific OGs were primarily salivary gland-specific in
early nodes but shifted to VG-specific at node 18 (Fig. 4c).
These shifts suggest that the organs of the digestive system
had to adapt to the new function, or loss of, the mid-
esophageal gland.

Discussion

Understanding how organs evolve is central to studying func-
tional innovation in animals. Our investigation into marine
gastropod mid-esophageal glands provides a unique perspec-
tive on the specific events driving the emergence of specialized
functions like venom production. Venom glands in animals

are highly specialized exocrine organs dedicated primarily
to toxin synthesis and secretion (Schendel et al. 2019).
Consistent with this, our functional enrichment analyses re-
veal that venom glands in cone snails are distinguished by ex-
ceptionally highly upregulation of protein secretion pathways,
consistent with findings in other lineages (Perry et al. 2020;
Barua and Mikheyev 2021; Zancolli et al. 2022). In contrast,
genes specifically expressed in the gland of Leiblein, and, to a
lesser extent, the esophageal gland, are associated with intra-
cellular digestion, particularly of proteins, as well as uptake
and storage of lipids and carbohydrates (Andrews and
Thorogood 2005). In particular, the esophageal gland appears
to play roles also in homeostasis, transport of fluids and nu-
trients, and secretion of digestive enzymes (Reid and Friesen
1980).

The phylogenetic distribution of these mid-esophageal
glands, combined with the functions inferred from their tran-
scriptomes, suggests a scenario of concerted functional and
morphological evolution closely linked to their trophic ecol-
ogy. Initially, the role of these glands in nutrient absorption
and intracellular digestion would have required prey to be
either predigested by salivary enzymes or naturally soft tissues. For example, species in the Muricidae family that drill holes into prey and consume soft tissues have a well-developed gland of Leiblein (Andrews and Thorogood 2005), and their salivary glands open near the tip of the proboscis, enabling immediate interaction of salivary products with prey (Kantor 1996). As gastropods diversified their diets and some species adopted macrophagous feeding strategies, the glands evolved accordingly, trending toward reduction or even loss. For instance, Buccinoidea species ingest prey whole or in large chunks (Kantor 2003) and possess a simple gland structure (Andrews and Thorogood 2005). In cases where the gland no longer functioned in absorption, it either disappeared, as in the Mitridae, or specialized for new functions, as seen in conoideans including cone snails. While it is often proposed that venom evolved to enable animals to capture and consume larger prey, our findings suggest an intriguing reversal in marine gastropods: An initial shift toward macrophagus feeding may have paved the way for venom evolution rather than vice versa. A comprehensive comparison of feeding strategies across lineages is needed to validate our proposed scenario.

The functional transformation of the mid-esophageal gland in Conoidea was facilitated by coordinated adaptations in other digestive organs, which shifted their functions accordingly. Our ancestral state reconstruction reveals extensive changes in gene expression across the phylogeny, reflecting the dynamic nature of this system and the adaptation of these snails to diverse ecological niches and feeding strategies. Notably, the venomous clade experienced more substantial expression changes and higher evolutionary rates than nonvenomous species, consistent with the drastic functional divergence of the venom gland from its homologous counterparts. However, the venom gland itself did not exhibit higher evolutionary rates than the homologous glands, nor higher rates compared to other organs. This suggests concerted evolution across the digestive systems, with other organs co-evolving to support the gland's new function. In the ancestor of the venomous species and the glandless *M. mitra* (node 18), several genes became downregulated in the gland while upregulated in the esophagus (Fig. 4a and c). Among these were genes involved in gluconeogenesis, such as galectin and PEPCK, indicating that the ancestral gland's role in producing energy from non-carbohydrate substrates was reassigned to other digestive tissues. This concerted adaptation is particularly evident in the glandless species, where many otherwise OEG- and LEG-specific genes are now esophagus-specific (Fig. 3b), especially those related to lysosomal functions. Mitridae snails feed on soft-bodied Sipuncula worms which they may ingest whole (Taylor 1989, 1993), or by pumping the worm's viscera (West 1990), including celomic fluids and eggs (West 1991), into the buccal cavity. In such cases, nutrients may be predigested by salivary enzymes, allowing the esophagus to take over functions previously performed by the gland of Leiblein.

During conoidean evolution, the mid-esophageal gland lost its original digestive functions but gained enhanced secretory capacity, primarily through the modulation of genes involved in pre-existing secretory pathways. Our analysis identified key upregulated genes in the ancestral venom gland (node 18) that encode proteins active in the ER. These include DAD1, which is critical for N-glycosylation and protein translocation (Zhang et al. 2016), disulfide isomerases, and a selenoprotein F likely involved in ER protein folding quality control (Takeda et al. 2014). Additionally, genes associated with the UPR and

ER stress response were enriched in cone snail venom glands, echoing patterns seen across other venomous lineages (Perry et al. 2020; Zancolli et al. 2022). This finding is significant for two reasons: First, gastropod venom glands reaffirm the trend of convergent transcriptomic evolution in venom glands across Metazoa (Zancolli et al. 2022). Second, the upregulation of these pathways appears unique to venom glands, as neither the extant homologous glands nor the ancestral organ shows this pattern based on transcriptome reconstruction. The strong upregulation of the UPR pathway in venom glands is notable. Even when compared to other exocrine organs with high secretory demands, such as salivary glands (this study) or the pancreas (Perry et al. 2020), venom glands show much higher expression levels. This indicates that heightened UPR expression is a distinct adaptation that evolved specifically within venom glands to support their unique physiological demands and, potentially, the direct control of toxin expression. Studies in snakes have led to a model where venom production activates the UPR, creating a positive feedback loop that enhances venom production through upregulation and binding of UPR transcription factors targeting toxin genes (Perry et al. 2020, 2022; Westfall et al. 2023).

In earlier nodes of the tree, OEG- and LEG-specific genes were expressed in the ancestral gland, while VG-specific genes were found in salivary glands—a pattern still evident in modern species (Figs. 3 and 4c). This convergence between venom and salivary glands is logical, as both are exocrine organs secreting products into the extracellular environment. Notably, the salivary glands of some neogastropods secrete toxins to immobilize prey (Ponte and Modica 2017), thereby expanding their function. In these cases, salivary secretions serve both endogenous roles (e.g. predigestion enzymes) and exogenous functions (e.g. toxins), indicating a functional convergence between these two nonhomologous organs. Interestingly, the salivary glands of the glandless species *M. mitra* share several tissue-specific genes with the venom glands. Unlike other snails, *Mitra* has salivary ducts that open at the tip of the epiproscis, an extendible muscular rod within the proboscis (Harasewych 2009), likely facilitating the direct delivery of secretions to prey (West 1990). We observed overexpression of cysteine-rich secreted proteins, peptidases, and serine proteases in *Mitra*'s salivary glands, which suggests roles in tissue degradation and potential toxin activity. Despite these functional similarities, salivary glands retained their original digestive role and only secondarily adopted a toxin-secreting function, while the venom gland became fully specialized solely for the latter. A plausible explanation for why the salivary glands did not evolve into specialized venom glands lies, again, in the diet—Mitridae snails primarily feed on Sipuncula worms, a relatively inactive prey group targeted by few other predators. This reduced competition may have lessened the selective pressure to develop a specialized venom apparatus, unlike cone snails that faced greater competition and adapted to prey on fast-moving organisms that required paralyzing toxins for successful capture.

## Conclusions

This study provides new insights into the genetic basis of functional innovation by examining the mid-esophageal glands of marine gastropods, with a particular focus on the evolution of venom production. Our results indicate that while mid-esophageal glands share a common origin, they have diverged in gene expression profiles and functions, shaped by
adaptations to different feeding strategies. In cone snails, an-
cestral digestive functions of the mid-esophageal gland were
relocated to other digestive tissues in a process of concerted
evolution, enabling the venom gland to specialize exclusively
in toxin production through modulation of pathways related
to secretion and cellular stress management. Overall, this
study underscores the link between transcriptome evolution
and functional divergence and identify the specific events lead-
ing to the emergence of new physiological functions.

**Materials and Methods**

Sample Collection and Sequencing

Between two and five individuals from Neogastropoda, and
two outgroup species, were collected in Koumac, New
Caledonia, under permit N°609011-55/2019/DEPART/JJC.
Individuals were kept for a minimum of 1 d and a maximum
of 3 d in aquaria with fresh sea water before dissection.
Salivary glands, foot, columellar muscle, esophagus, esopha-
geal gland, gland of Leiblein, dorsal glandular folds (only in
*M. tenuirostrum*), venom gland, and muscular bulb (only in
cone snails) were dissected and preserved in *RNAlater*
(Invitrogen) (supplementary dataset S1, Supplementary
Material online).

Tissue samples were homogenized with Trizol (Invitrogen) and total RNA purified with the PureLink RNA Mini kit (Invitrogen) with an additional DNase I treatment following manufacturer's protocol. cDNA libraries were constructed using the NEBnext Ultra II Directional RNA Library Preparation Kit with polyA selection and dUTP method (New England BioLabs) and sequenced on an Illumina NovaSeq 2×150 bp. Raw paired-end reads were checked with FastQC 0.11.9 (Andrews 2010), quality-filtered and trimmed with FastP 0.20.1 (Chen et al. 2018).

De Novo Assembly and Annotation

All the reads from a species were pooled to generate de novo
transcriptome assemblies using rnaSpades 3.15.2
(Bushmanova et al. 2019). To reduce assembly size and redun-
dancy, and remove spurious transcripts, we adopted a series of
filtering steps. First, we translated the transcripts to amino
acid sequences using Borf 1.2.1 (Signal and Kahlke 2021)
and kept only those with a complete open read frame.
Second, we compared the translated sequences with BlastP
(Altschul et al. 1997) against a suite of databases (downloaded
on 30.07.2021) including UniprotKB/Swiss-Prot (The UniProt
Consortium 2023), a set of 11 Gastropoda genomes
(supplementary table S2, Supplementary Material online),
Conoserver (Kaas et al. 2012), and Tox-Prot (Jungo et al.
2012). Protein domains were identified with PSIblast
(Altschul et al. 1997) against the Pfam (Mistry et al. 2021)
and Cdd (Wang et al. 2023) databases. Only hits with eval-
ue < 1e-05 against at least one database were retained.
Additionally, we annotated signal peptides with signalP 6.0
(Teufel et al. 2022). Putative conotoxins were predicted with
the ConoPrec tool available on ConoServer (Kaas et al.
2012). All sequences with the UniprotKB/Swiss-Prot top hit
(i.e. lowest evalue) to a non-Metazoa organism were removed.
We then trimmed the retained transcripts to their coding re-
gion, reduced redundancy by clustering sequences with 98%
or more identical nucleotide sequences with CD-HIT-EST
4.6 (Li and Godzik 2006), and kept only transcripts with
more than one read count in at least one library. Finally, we

evaluated assembly completeness using Omark on the web-
server omark.omabrowser.org (Nevers et al. 2025).

For each species, we performed GO annotation by combin-
ing the annotation from Pannzer2 (Törönen et al. 2018) and
DeepGOPlus 1.0.2 (Kulmanov and Hoehndorf 2020). Based
on the score values distribution, we used 0.3 as a score thresh-
old for both methods.

Within-Species Analysis

Expression Levels

We mapped all libraries to the respective species assembly with
Kallisto 0.48.0 (Bray et al. 2016) with 100 bootstrap and
quantified gene expression using the package Sleuth 0.30.1
(Pimentel et al. 2017) using R 4.2.2 (R Core Team 2019).
We employed a suite of quality control steps to identify and re-
move outliers. First, we analyzed read count distribution with
*vioplot* 0.4.0 and removed samples with particularly different
distributions, generally with read counts lower than the aver-
age. Then, we utilize dimensionality reduction techniques, in-
cluding PCA (*dudi_PCA*) and multidimensional scaling
(*plotMDS*) to assess whether samples were clustering by tissue
type. Normalized estimated counts and TPM abundances
were then re-calculated after outlier libraries’ removal. As in
some species samples still tended to cluster by individual ra-
ther than tissue type, we account for the specimen effect using
an empirical Bayes method implemented via the *ComBat_seq*
function in the package sva 3.46.0 (Leek et al. 2012) which
specifically targets RNA-seq data. *ComBat_seq* uses a nega-
tive binomial regression model that retains the integer nature
of count data of RNA-seq, resulting in better statistical power
and control of false positives compared to data adjusted by
other methods (Zhang et al. 2020). We run *ComBat_seq*
with the following parameters: *batch* = specimen, *group* = tis-
*sue_type, full_mod* = TRUE.

Tissue-Specific Gene Sets

For each species, we identified tissue-specific genes based on
their fold change (FC) calculated as the ratio between the
TPM value of the first most highly expressed tissue and the
TPM value of the second most highly expressed tissue. A
gene with TPM≥2 and FC≥2 was classified as specific of
the top tissue. We validated our FC method by confirming
that the tissue-specific genes identified were also significantly
differentially expressed when analyzed using the likelihood ra-
tio test in Sleuth. Differences in the number of tissue-specific
genes between glands were tested by means of one-tailed
t-tests. Functional enrichment of tissue-specific gene sets was
performed with TopGO 2.50.0 (Alexa and Rahnenfuhrer
2019) using the *elim* algorithm and Fisher test. The fore-
ground was the list of tissue-specific genes while the back-
ground included all the genes expressed in that species.

Between-Species Analysis

Orthogroup Expression Matrix

Amino acid sequences were assigned to orthogroups (OGs)
with the OrthoDB standalone pipeline OrthoLoger 3.0.2
(Kuznetsov et al. 2023) using default parameters. Since most
OGs included more than one gene per species (i.e.
one-to-many or many-to-many orthologs), we randomly se-
lected one representative sequence for each OG in each species
and used the TPM value estimated for that gene as the
orthogroup expression value. This method was shown to be
robust (Zancolli et al. 2022). Alternatively, we calculated
the mean TPM values across all genes belonging to the same
OG. All samples were then merged into a multispecies multi-
tissue matrix and the expression levels corrected for the species
effect using *ComBat_seq* (Zhang et al. 2020). All the down-
stream analyses reported in the main text are based on the ran-
dom matrix, while the results from the mean-based matrix are
reported in Supplementary Material online.

To have an overview of global transcriptome similarity
across tissues and species, we calculated pairwise distances
as 1- Spearman correlation and used it to reconstruct a gene
expression tree using the neighbor-joining method.
Additionally, we calculated the proportion of shared OGs
among tissue-specific genes. For downstream analyses, OGs
were classified tissue-specific if at least one gene within that
OG was specific in a tissue, and if it was found specific in at
least two species.

Gene Expression Evolution Analysis

We analyzed changes in gene expression with the program CAGEE 1.1.1 (Bertram et al. 2023) which uses Brownian motion to model gene expression across a phylogenetic tree. The tree was derived from the phylogeny of Fedosov et al. (2024) by excluding the families not encompassed in this study and retaining only the branches most closely related to the species that we examined. The rates of expression changes ($\sigma^2$) were calculated for the mid-esophageal glands, salivary glands, and esophagus separately. We fit a series of nested models in which $\sigma^2$ varies across branches of the species tree to test for different hypotheses as outlined in the results. Ancestral transcriptomes at inner nodes reconstructed in the best fit model were used to calculate ancestral tissue specificity for each OG as we did for extant species.

**Supplementary Material**

*Supplementary material* is available at *Molecular Biology and Evolution* online.

Acknowledgments

We are grateful to Philippe Buchet, the team of the expeditions
Koumac 2.1 and 2.3 of the Muséum National d'Histoire
Naturelle, Paris, and the Province Nord and Government of
New Caledonia for providing the logistic support and expert-
ise; David Massemin and Patrick Marti for their skills in find-
ing the species needed. We thank Alexander Fedesov for
providing the phylogenetic tree, Consolée Aletti for support
with the RNA extraction, the personnel of the Genomic
Technologies Facility and Genewiz for the RNA-seq experi-
ment, and members of the Robinson Rechavi's lab for bio-
informatics support and scientific discussions. The first
version of the manuscript was greatly improved thanks to
the feedback of two anonymous reviewers and six additional
reviewers through Qeios.

Author Contributions

G.Z., M.V.M., N.P., Y.K., and M.R.-R. conceived the study
design; Y.K. and M.V.M. perform the dissections; G.Z. gener-
ated the RNA-seq data, performed the data analysis, and
wrote the manuscript; A.B. performed the GO annotation
and G.C. the orthogroup assignment. All the authors

contributed to the interpretation and discussion of the results,
and to the final version of the manuscript.

Funding

This work has received funding from the European Union's Horizon 2020 research and innovation programme under Marie Skłodowska-Curie grant agreement 845674 to G.Z., and grant agreement 865101 to N.P.

Conflict of Interest

The authors declare no competing interests.

Data Availability

The RNA-seq data generated in this study have been deposited
in the NCBI SRA archive with the accession number
PRJNA1158673. Additional data generated in this study are
available at 10.5281/zenodo.13685166. Additional informa-
tion are provided as Supplementary files.

References

Abalde S, Tenorio MJ, Afonso CML, Zardoya R. Conotoxin diversity in *Chelyconus ermineus* (Born, 1778) and the convergent origin of piscivory in the Atlantic and Indo-Pacific cones. *Genome Biol Evol.* 2018;10(10):2643–2662. https://doi.org/10.1093/gbe/evy150.

Abdelkrim J, Aznar-Cormano L, Fedosov AE, Kantor YI, Lozouet P, Phuong MA, Zaharias P, Puillandre N. Exxon-capture-based phylogeny and diversification of the venomous gastropods (Neogastropoda, Conoidea). *Mol Biol Evol.* 2018;35(10): 2355–2374. https://doi.org/10.1093/molbev/msy144.

Alexa A, Rahnenfuhrer J. topGO: Enrichment analysis for gene ontol-
ogy. R package version 2.38.1. 2019.

Altschul SF, Madden TL, Schäffer AA, Zhang J, Zhang Z, Miller W, Lipman DJ. Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. *Nucleic Acids Res.* 1997;25(17):3389–3402. https://doi.org/10.1093/nar/25.17.3389.

Andrews EB, Thorogood KE. An ultrastructural study of the gland of Leiblein of muricid and nassariid neogastropods in relation to function, with a discussion on its homologies in other caenogastropods. *J Molluscan Stud.* 2005;71(3):269–300. https://doi.org/10.1093/mollus/eyi036.

Andrews S. FastQC: a quality control tool for high throughput sequence
data. 2010. https://www.bioinformatics.babraham.ac.uk/projects/
fastqc/

Barua A, Mikheyev AS. An ancient, conserved gene regulatory network
led to the rise of oral venom systems. *Proc Natl Acad Sci U S A.*
2021;118(14):e202131118. https://doi.org/10.1073/pnas.202131118

Benton MJ, Dhouailly D, Jiang B, McNamara M. The early origin of feathers. *Trends Ecol Evol.* 2019;34(9):856–869. https://doi.org/10.1016/j.tree.2019.04.018.

Bertram J, Fulton B, Tourigny JP, Peña-Garcia Y, Moyle LC, Hahn MW. CAGEE: computational analysis of gene expression evolution. *Mol Biol Evol.* 2023;40(5):msad106. https://doi.org/10.1093/molbev/msad106.

Bray NL, Pimentel H, Melsted P, Pachter L. Near-optimal probabilistic
RNA-seq quantification. *Nat Biotechnol.* 2016;34(5):525–527.
https://doi.org/10.1038/nbt.3519.

Bushmanova E, Antipov D, Lapidus A, Prjibelski AD. rnaSPAdes: a *de novo* transcriptome assembler and its application to RNA-Seq data. *GigaScience.* 2019;8(9):giz100. https://doi.org/10.1093/gigascience/giz100.
Casewell NR, Wüster W, Vonk FJ, Harrison RA, Fry BG. Complex cocktails: the evolutionary novelty of venoms. *Trends Ecol Evol.* 2013;28(4):219–229. https://doi.org/10.1016/j.tree.2012.10.020.

Chen S, Zhou Y, Chen Y, Gu J. fastp: an ultra-fast all-in-one FASTQ preprocessor. *Bioinformatics* (Oxford, Engl). 2018;34(17):i884–i890. https://doi.org/10.1093/bioinformatics/bty560.

Dutt M, Dutertre S, Jin A-H, Lavergne V, Alewood PF, Lewis RJ. Venomics reveals venom complexity of the piscivorous cone snail, *Conus tulipa*. *Mar Drugs*. 2019;17(1):71. https://doi.org/10.3390/ md17010071.

Fedosov A, Tucci CF, Kantor Y, Farhat S, Puillandre N. Collaborative expression: transcriptomics of *Conus virgo* suggests contribution of multiple secretory glands to venom production. *J Mol Evol*. 2023;91(6):837–853. https://doi.org/10.1007/s00239- 023-10139-8.

Fedosov AE, Zaharias P, Lemarcis T, Modica MV, Holford M, Oliverio M, Kantor YI, Puillandre N. Phylogenomics of Neogastropoda: the backbone hidden in the bush. *Syst Biol*. 2024;73(3):521–531. https://doi.org/10.1093/sysbio/syae010.

Harasewych MG. Anatomy and biology of *Mitra cornea* Lamarck, 1811 (Mollusca, Caenogastropoda, Mitridae) from the Azores. *Açoreana*. 2009;6:121–135.

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(1):15–25. https://doi.org/10.1093/molbev/msn211.

Jungo F, Bougueleret L, Xenarios I, Poux S. The UniProtKB/Swiss-Prot Tox-Prot program: a central hub of integrated venom protein data. *Toxicon*. 2012;60(4):551–557. https://doi.org/10.1016/j.toxicon. 2012.03.010.

Kaas Q, Yu R, Jin A-H, Dutertre S, Craik DJ. ConoServer: updated content, knowledge, and discovery tools in the conopeptide database. *Nucleic Acids Res.* 2012;40(D1):D325–D330. https://doi.org/10. 1093/nar/gkr886.

Kantor Y. Phylogeny and relationships of Neogastropoda. In: Taylor J, editor. *Origin and evolutionary radiation of the Mollusca*. Oxford: Oxford University Press; 1996. p. 221–230.

Kantor YI. Comparative anatomy of the stomach of Buccinoidea (Neogastropoda). *J Molluscan Stud*. 2003;69(3):203–220. https:// doi.org/10.1093/mollus/69.3.203.

Kulmanov M, Hoehndorf R. DeepGOLive: improved protein function prediction from sequence. *Bioinformatics*. 2020;36(2):422–429. https://doi.org/10.1093/bioinformatics/btz595.

Kuznetsov D, Tegenfeldt F, Manni M, Seppey M, Berkeley M, Kriventseva EV, Zdobnov EM. OrthoDB v11: annotation of orthologs in the widest sampling of organismal diversity. *Nucleic Acids Res.* 2023;51(D1):D445–D451. https://doi.org/10.1093/nar/ gkac998.

Leek JT, Johnson WE, Parker HS, Jaffe AE, Storey JD. The sva package for removing batch effects and other unwanted variation in high-throughput experiments. *Bioinformatics*. 2012;28(6):882–883. https://doi.org/10.1093/bioinformatics/bts034.

Li W, Godzik A. CD-HIT: a fast program for clustering and comparing large sets of protein or nucleotide sequences. *Bioinformatics*. 2006;22(13):1658–1659. https://doi.org/10.1093/bioinformatics/ btl158.

Lobo-da-Cunha A. Structure and function of the digestive system in molluscs. *Cell Tissue Res.* 2019;377(3):475–503. https://doi.org/ 10.1007/s00441-019-03085-9.

Mistry J, Chuguransky S, Williams L, Qureshi M, Salazar GA, Sonnhammer ELL, Tosatto SCE, Paladin L, Raj S, Richardson LJ, et al. Pfam: the protein families database in 2021. *Nucleic Acids Res.* 2021;49(D1):D412–D419. https://doi.org/10.1093/nar/ gkaa913.

Modica MV, Holford M. 2010. The Neogastropoda: evolutionary innovations of predatory marine snails with remarkable pharmacological potential. In: Pontarotti P, editor. *Evolutionary biology—concepts, molecular and morphological evolution: 13th Meeting 2009*. Berlin, Heidelberg: Springer. p. 249–270. https://doi.org/10.1007/ 978-3-642-12340-5_15.

Morales Duque H, Campos Dias S, Franco OL. Structural and functional analyses of cone snail toxins. *Mar Drugs*. 2019;17(6):370. https:// doi.org/10.3390/md17060370.

Nevers Y, Warwick Vesztyrocy A, Rossier V, Train C-M, Altenhoff A, Dessimoz C, Glover NM. Quality assessment of gene repertoire annotations with OMArk. *Nat Biotechnol.* 2025;43(1):124–133. https://doi.org/10.1038/s41587-024-02147-w.

Nguyen LTT, Craik DJ, Kaas Q. Bibliometric review of the literature on cone snail peptide toxins from 2000 to 2022. *Mar Drugs*. 2023;21(3):154. https://doi.org/10.3390/ md21030154.

Olivera BM, Watkins M, Bandyopadhyay P, Imperial JS, de la Cotera EPH, Aguilar MB, Vera EL, Concepcion GP, Lluisma A. Adaptive radiation of venomous marine snail lineages and the accelerated evolution of venom peptide genes. *Ann N Y Acad Sci.* 2012;1267(1):61–70. https://doi.org/10.1111/j.1749-6632.2012.06603.x.

Page LR. Development and evolution of adult feeding structures in Caenogastropods: overcoming larval functional constraints. *Evol Dev.* 2000;2(1):25–34. https://doi.org/10.1046/j.1525-142x.2000. 00017.x.

Page LR. Larval and metamorphic development of the foregut and proboscis in the caenogastropod *Marsenina (Lamellaria) stearnsii*. *J Morphol.* 2002;252(2):202–217. https://doi.org/10.1002/jmor. 1099.

Page LR. Metamorphic remodeling of a planktotrophic larva to produce the predatory feeding system of a cone snail (Mollusca, Neogastropoda). *Biol Bull.* 2011;221(2):176–188. https://doi.org/ 10.1086/BBLv221n2p176.

Page LR. Developmental modularity and phenotypic novelty within a biphasic life cycle: morphogenesis of a cone snail venom gland. *Proc Biol Sci.* 2012;279(1726):77–83. https://doi.org/10.1098/ rspb.2011.0501.

Page LR, Hookham B. The gastropod foregut—evolution viewed through a developmental lens. *Can J Zool.* 2017;95(4):227–238. https://doi.org/10.1139/cjz-2016-0194.

Perry BW, Gopalan SS, Pasquesi GIM, Schield DR, Westfall AK, Smith CF, Koludarov I, Chippindale PT, Pellegrino MW, Chuong EB, et al. Snake venom gene expression is coordinated by novel regulatory architecture and the integration of multiple co-opted vertebrate pathways. *Genome Res.* 2022;32(6):1058–1073. https://doi.org/ 10.1101/gr.276251.121.

Perry BW, Schield DR, Westfall AK, Mackessy SP, Castoe TA. Physiological demands and signaling associated with snake venom production and storage illustrated by transcriptional analyses of venom glands. *Sci Rep.* 2020;10(1):18083. https://doi.org/10. 1038/s41598-020-75048-y.

Pimentel H, Bray NL, Puente S, Melsted P, Pachter L. Differential analysis of RNA-seq incorporating quantification uncertainty. *Nat Methods.* 2017;14(7):687–690. https://doi.org/10.1038/nmeth. 4324.

Ponte G, Modica MV. Salivary glands in predatory mollusks: evolutionary considerations. *Front Physiol.* 2017;8:580. https://doi.org/10. 3389/fphysiol.2017.00580.

R Core Team. *R: a language and environment for statistical computing*. Vienna: R Foundation for Statistical Computing; 2019. http://www. R-project.org/

Reid R, Friesen J. The digestive system of the moon snail *Polinices lewisi* (Gould, 1847) with emphasis on the role of the oesophageal gland. *Veliger*. 1980;23:25–34.

Safavi-Hemami H, Young ND, Williamson NA, Purcell AW. Proteomic interrogation of venom delivery in marine cone snails: novel insights into the role of the venom bulb. *J Proteome Res.* 2010;9(11): 5610–5619. https://doi.org/10.1021/pr100431x.

Salisbury SM, Martin GG, Kier WM, Schulz JR. Venom kinematics during prey capture in *Conus*: the biomechanics of a rapid injection system. *J Exp Biol.* 2010;213(5):673–682. https://doi.org/10.1242/jeb. 035550.

Schendel V, Rash LD, Jenner RA, Undheim EAB. The diversity of venom: the importance of behavior and venom system morphology in understanding its ecology and evolution. *Toxins (Basel).* 2019;11(11):666. https://doi.org/10.3390/toxins11110666.
Shiomi K, Kawashima Y, Mizukami M, Nagashima Y. Properties of proteinaceous toxins in the salivary gland of the marine gastropod (*Monoplex echo*). *Toxicon*. 2002;40(5):563-571. https://doi.org/10.1016/S0041-0101(01)00256-2.

Signal B, Kahlke T. 2021. Borf: improved ORF prediction in de-novo assembled transcriptome annotation. *bioRXiv* 439551. https://doi.org/10.1101/2021.04.12.439551, preprint: not peer reviewed.

Takeda Y, Seko A, Hachisu M, Daikoku S, Izumi M, Koizumi A, Fujikawa K, Kajihara Y, Ito Y. Both isoforms of human UDP-glucose:glycoprotein glucosyltransferase are enzymatically active. *Glycobiology*. 2014;24(4):344-350. https://doi.org/10.1093/glycob/cwt163.

Taylor JD. The diet of coral-reef Mitridae (Gastropoda) from Guam; with a review of other species of the family. *J Nat Hist*. 1989;23(1): 261-278. https://doi.org/10.1080/00222938900770141.

Taylor JD. 1993. Dietary and anatomical specialization of mitrid gastropods (Mitridae) at Rottnest island, Western Australia. In: Wells FE, Walker DI, Kirkman H, Lethdridge R, editors. *Proceedings of the fifth international marine biological workshop: the marine flora and fauna of Rottnest Island, Western Australia*. Perth: Western Australian Museum. p. 583-599.

Teufel F, Almagro Armenteros JJ, Johansen AR, Gíslason MH, Pihl SI, Tsirigos KD, Winther O, Brunak S, von Heijne G, Nielsen H. SignalP 6.0 predicts all five types of signal peptides using protein language models. *Nat Biotechnol*. 2022;40(7):1023-1025. https://doi.org/10.1038/s41587-021-01156-3.

The UniProt Consortium. UniProt: the universal protein knowledgebase in 2023. *Nucleic Acids Res*. 2023;51(D1):D523-D531. https://doi.org/10.1093/nar/gkac1052.

Törönen P, Medlar A, Holm L. PANNZER2: a rapid functional annotation web server. *Nucleic Acids Res*. 2018;46(W1):W84-W88. https://doi.org/10.1093/nar/gky350.

Wang J, Chitsaz F, Derbyshire MK, Gonzales NR, Gwadz M, Lu S, Marchler GH, Song JS, Thanki N, Yamashita RA, et al. The conserved domain database in 2023. *Nucleic Acids Res*. 2023;51(D1):D384-D388. https://doi.org/10.1093/nar/gkac1096.

West TL. Feeding behavior and functional morphology of the epiproscos of *Mitra idaea* (Mollusca: Gastropoda: Mitridae). *Bull Mar Sci*. 1990;46:761-779.

West TL. Functional morphology of the proboscis of *Mitra catalinae* Dall 1920 (Mollusca: Gastropoda: Mitridae), and the evolution of the mitrid epiproscosis. *Bull Mar Sci*. 1991;48:702-718.

Westfall AK, Gopalan SS, Perry BW, Adams RH, Saviola AJ, Mackessy SP, Castoe TA. Single-cell heterogeneity in snake venom expression is hardwired by co-option of regulators from progressively activated pathways. *Genome Biol Evol*. 2023;15(6):evad109. https://doi.org/10.1093/gbe/evad109.

Zancolli G, Casewell NR. Venom systems as models for studying the origin and regulation of evolutionary novelties. *Mol Biol Evol*. 2020;37(10):2777-2790. https://doi.org/10.1093/molbev/msaa133.

Zancolli G, Reijnders M, Waterhouse RM, Robinson-Rechavi M. Convergent evolution of venom gland transcriptomes across Metazoa. *Proc Natl Acad Sci U S A*. 2022;119(1):e2111392119. https://doi.org/10.1073/pnas.2111392119.

Zhang Y, Cui C, Lai Z-C. The defender against apoptotic cell death 1 gene is required for tissue growth and efficient N-glycosylation in *Drosophila melanogaster*. *Dev Biol*. 2016;420(1):186-195. https://doi.org/10.1016/j.ydbio.2016.09.021.

Zhang Y, Parmigiani G, Johnson WE. ComBat-seq: batch effect adjustment for RNA-seq count data. *NAR Genom Bioinform*. 2020;2(3):lqaa078. https://doi.org/10.1093/nargab/lqaa078.

Zuniga A. Next generation limb development and evolution: old questions, new perspectives. *Development*. 2015;142(22):3810-3820. https://doi.org/10.1242/dev.125757.