Three adult specimens of C. eburneus (Ce1, Ce2, and Ce3) and C. tessulatus (Ct1, Ct2, and Ct3) were collected by the contracted local fishermen in Caw-oy Lapu-Lapu City, Cebu, Philippines. The specimens were initially identified by morphological examination of the shell. Prior to dissection, live snails were acclimatized for 24 hours in an aerated (improvised) aquarium filled with sea water. The snails were dissected directly on ice, and venom ducts were carefully obtained and stored in 1 ml of RNAlater^® (Invitrogen) at 4°C prior to long term storage at -80°C. The total RNA was isolated using the Trizol extraction method. Briefly, thawed venom duct was homogenized using 2.0 mm of ZR bashing beads (ZYMO research) in a bead beater (Precellys, Berlin Technologies) with 1 ml TRIzol reagent (Sigma-Aldrich TRI reagent^®). The total RNA was extracted from the homogenate using chloroform and precipitated using isopropanol. The pellet was washed twice with 75% EtOH and air dried. The pelleted RNA was resuspended in 50 ul of RNAase-free water, and clean-up was done using the Qiagen purification kit following the manufacturers protocol. The quality and quantity of total RNA was evaluated using the Agilent TapeStation and Qubit, respectively ( Supplementary Figure S5 ). Individual library was constructed in each sample using the TruSeq Stranded mRNA Library prep kit and the sequencing was performed in Illumina NextSeq550 using Mid Output v2.5 with 300 cycles. All library preparations and sequencing were done in the sequencing facility of PCARI-SGCL, University of the Philippines.

Six sets of paired-end sequence reads from three individuals of C. eburneus and C. tessulatus were pre-processed mainly using the suggested best practices for de novo transcriptome assembly with Trinity by the Harvard FAS Informatics group (https://informatics.fas.harvard.edu/best-practices-for-de-novo-transcriptome-assembly-with-trinity.html). Briefly, the initial and final quality assessment of raw and filtered sequence data were done using Fastqc (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Erroneous base calls were then corrected based on k-mer composition using Rcorrector (Song and Florea, 2015), and reads that were deemed unfixable were subsequently removed using the FilterUncorrectabledPEfastq.py script obtained from the Transcriptome Assembly Tools repository of the Harvard Informatics group (https://github.com/harvardinformatics/TranscriptomeAssemblyTools). Adapter and quality trimming of the error corrected reads was implemented afterwards using the tool Trim Galore (https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/). The remaining reads were further mapped against the SILVA database of rRNA sequences (Quast et al., 2013) using Bowtie2 (Langmead and Salzberg, 2012), and only reads that did not map to any of the rRNA sequences were retained. For each of the individual samples, >90% of the reads passed the filtering step ( Supplementary Table S7 ).

Quality filtered paired-end read pairs from three individuals of each species were co-assembled using three different tools: Trinity (Haas et al., 2013) implemented with k-mers 25 and 31, Trans-ABySS (Robertson et al., 2010), and rnaSPAdes (Bushmanova et al., 2019) – the latter two using default parameters. The resulting assemblies from each of the tools and parameters were then combined and clustered at 100% sequence similarity threshold using CD-HIT (Fu et al., 2012).

Transcripts from the resulting combined assemblies (one each for C. eburneus and C. tessulatus) were further assessed and filtered using TransRate (Smith-Unna et al., 2016), and only those included in the set of good contigs were used in the subsequent analyses. From the set of good contigs, transcript expression levels for each of the individual samples were estimated using Salmon (Patro et al., 2017).

Open reading frames (ORFs) were detected from each of the assembled transcripts using Transdecoder (https://github.com/TransDecoder/TransDecoder), which also outputs the corresponding amino acid sequences for the predicted coding regions. Note that only the best coding region per transcript was retained, either by being the longest predicted ORF for that particular transcript or by having significant sequence similarities with known conopeptides as determined by running the blastp algorithm of the Basic Local Alignment and Search Tool (BLAST) (Altschul et al., 1990) against a custom curated database of conopeptide sequences collected from the following data sources: GenBank (National Center for Biotechnology Information, USA), UniProt (Swiss Institute of Bioinformatics, Switzerland), and ConoServer (Kaas et al., 2012).

Amino acid sequences of transcripts with conopeptide matches were further analyzed using the ConoPrec tool in ConoServer (http://www.conoserver.org/?page=conoprec) in order to determine the precursor peptide structure, particularly the signal and mature peptide sequences. The same tool also assigns gene superfamily and cysteine framework classifications to each of the putative conopeptide sequences.

Putative conopeptides with known gene superfamily or cysteine framework classifications were gathered together with their corresponding expression estimates. Scatter plots of cumulative expression levels per gene superfamily or cysteine framework plotted against the number of unique conopeptide sequences for each of the given classifications were generated using custom scripts implemented in R (R Core Team, 2021). Venn diagrams of shared conopeptides within and between species were generated also using R.

In order to identify possibly novel gene superfamilies, signal sequences for each of the putative conopeptides as predicted by ConoPrec were clustered via agglomerative clustering at 75% similarity threshold with minimum distance single linkage using USEARCH (Edgar, 2010). The resulting cluster tree was then viewed and manually inspected using FigTree (http://tree.bio.ed.ac.uk/software/figtree/). Cluster of signal sequences with unknown gene superfamily having at least four member sequences were then classified as putative novel superfamilies. For each putative novel superfamily cluster, member signal sequences were aligned using the L-INS-I algorithm in MAFFT (Katoh et al., 2002). The resulting alignments were then used to generate sequence profiles using WebLogo (Crooks et al., 2004).

Putative conopeptide sequences without known cysteine framework classifications were analyzed for sequence similarity with previously identified disulfide-poor conopeptide sequences in the ConoServer database using the blastp algorithm in BLAST. In this case, conopeptides are deemed disulfide-poor if they belong to any of the following classes: conopressin, contryphan, conoCAP, conoGAY, conantokin, contulakin, conorfamide, conophan, conomap, conomarphin, conolysin, cono-NPY, and hormone-like including insulin (Lebbe and Tytgat, 2016). Alignment of putative conantokin precursor sequences classified under the B1 conopeptide superfamily was also done and viewed using Aliview (Larsson, 2014).

Predicted mature peptide sequences were used as BLAST queries against a database of mature conopeptide sequences obtained from Conus species with known diet preference based on ConoServer information. Sequences from C. eburneus and C. tessulatus were excluded from the custom database in order to get the nearest conopeptide match from a different species.

Sequences of delta conopeptides, believed to be involved in the transition from vermivorous to piscivorous diet, were also clustered using USEARCH (75% similarity threshold, minimum distance single linkage). The delta conopeptides from this study, as well as those from the UniProt database, with precursor sequence lengths of at least 40 a.a. were included in the analysis. Agglomerative clustering was done for both full precursor and predicted mature sequences, and the resulting cluster trees were viewed using FigTree.

From the non-redundant merged assemblies, a total of 201,656 good quality contigs were obtained from the transcriptome of C. eburneus and 282,501 from C. tessulatus. Among these assembled transcripts, 82,170 were identified to contain candidate coding regions from C. eburneus and 119,373 from C. tessulatus; of which, only 33,414 and 51,981 were predicted to have complete open reading frames, respectively. Nevertheless, based on the statistics reported by the transcript abundance estimation tool Salmon, an overall mapping rate of around 97% was observed for the quasi-mapping of sequence reads from all individual samples against the assemblies corresponding to their respective species. The pertinent assembly metrics are shown in Table 1 .

Sequence similarity search against a custom curated database, composed of conopeptide sequences from NCBI GenBank, UniProt and ConoServer, identified significant matches with 1,517 C. eburneus transcripts and 1,640 C. tessulatus transcripts, which corresponds to 954 and 1,063 unique peptide sequences, respectively. Among these, 424 sequences from C. eburneus were tagged to have complete coding regions, with 260 of these sequences predicted to carry a known conopeptide superfamily or cysteine framework. For C. tessulatus, 537 were tagged to have complete coding regions, with known superfamilies or cysteine frameworks identified in a subset of 339 sequences.

Note that a preliminary set of proteomic data from C. eburneus and C. tessulatus (unpublished) yielding a set of 26 and 20 peptide fragments, respectively, were all mapped to 22 unique transcript sequences from this study (8 from C. eburneus and 14 from C. tessulatus). However, because only a small fraction of the conopeptide sequences can be observed in the proteome due to the sheer throughput of transcriptomic analysis, we would like to caution that some of the transcripts reported in this study might not be translated into functional venom peptides. Similar to most other transcriptomic studies, the presence of assembly artifacts cannot be fully discounted as well. Nevertheless, we only included high confidence precursor conopeptide predictions in subsequent downstream analyses by filtering for those with significant sequence similarities with known conopeptides, as well as those with either a known gene superfamily or cysteine framework.

Among those with predicted complete coding sequences, the set of putative conopeptide transcripts with known superfamilies or cysteine frameworks were observed to be mostly different between C. eburneus and C. tessulatus. Based on sequence similarities of the predicted mature peptide regions, only 23 sequences overlap between the two sister species at 70% identity threshold, 12 of which have matching regions that are 100% identical ( Supplementary Figure S1 ). The comparison of estimated expression levels, predicted cysteine frameworks and gene superfamilies, as well as the actual peptide sequences of these 12 highly identical conopeptides shared by both Tesseliconus species are shown in Table 2 . Six of these putative conopeptides were predicted to have unknown, possibly novel, gene superfamily classifications.

Using the number of unique precursor conopeptide sequences as diversity estimate, it appears that conopeptides with cysteine framework VI/VII had the highest diversity in both C. eburneus and C. tessulatus individuals sampled in this study ( Figures 1A, B ). However, in terms of cumulative average expression, framework XIV was observed to be more abundantly expressed than VI/VII in both species. The majority of conopeptides predicted to have the XIV cysteine framework have unknown superfamily classifications – some of which have signal sequences forming exclusive clusters of possibly novel gene superfamilies (Clusters 3, 7 and 15), as described further in the latter sections of this paper, whereas those with superfamily classifications all fall under the A superfamily. Notably, the most highly expressed cysteine framework in C. eburneus, framework III, had substantially lower expression and diversity in C. tessulatus.

In terms of gene superfamily classification, O1 had the highest diversity in C. eburneus and second highest in C. tessulatus ( Figures 1C, D ). For C. tessulatus, the highest diversity and expression was observed for the O2 superfamily, whereas the most abundantly expressed superfamily in C. eburneus was M. Other notable differences include those of the A and T superfamilies whose expression levels were substantially higher in C. tessulatus than in C. eburneus, as well as the P superfamily whose representatives were detected exclusively in C. tessulatus.

Among the 260 predicted conopeptides with known cysteine framework or gene superfamily in C. eburneus, we found that 102 (39.23%) were expressed in all three individuals of the said species included in this study. An additional 91 (35.00%) were shared by two individuals, and 67 (25.77%) were uniquely expressed ( Figure 2A ). A relatively similar pattern was observed for C. tessulatus, with 142 (41.89%) of the 339 predicted conopeptides with known superfamily or framework found to be expressed in all three individuals, 110 (32.45%) shared by two samples, and 87 (25.66%) were unique to a particular individual ( Figure 2B ).

Although most of the conopeptides detected were shared by individuals of the same species, their relative expression levels appear to differ. Supplementary Figures S2 , S3 show the cumulative expression of conopeptides per cysteine framework or gene superfamily for both C. eburneus and C. tessulatus, respectively. For C. eburneus, individual 3 (Ce3) registered a relatively different pattern of cumulative expression, particularly for peptides with cysteine frameworks VI/VII (higher expression) and XIV (lower expression), compared to those observed in Ce1 and Ce2. Similarly, C. tessulatus individual Ct2 was observed to have higher VI/VII and lower XIV cumulative expression relative to Ct1 and Ct3. Other notable differences in expression levels across individuals of the same species include the following: O2 and T superfamilies and V cysteine framework in C. eburneus, as well as the T and M superfamilies and the III, V and XVI cysteine frameworks in C. tessulatus.

Agglomerative clustering (75% similarity threshold) of predicted conopeptide signal sequences from both C. eburneus and C. tessulatus revealed 10 known superfamily clusters and 18 possibly novel gene superfamilies with at least four representative Tesseliconus conopeptide sequences ( Figure 3 ). The O1 and M superfamilies form the largest sequence clusters, while O2, T, A and O3 also have relatively large representations. On the other hand, clusters 1, 3, 4, and 6 contain the highest number of protein precursors among the novel clusters. Table 3 lists the designated major gene superfamily clusters (both known and novel), the cysteine frameworks associated with these superfamilies, and the number of precursor peptides comprising these clusters. Sequences comprising the potentially novel gene superfamily clusters are also indicated in Supplementary Tables 1 , 2 , designated as cluster_01 to cluster_18 indicating their tentative numeric superfamily classifications.

Relatively larger superfamily clusters with at least 20 precursor peptides were mainly associated with multiple (three to four) cysteine frameworks, except for the O1 superfamily that has been associated exclusively with framework VI/VII. We do note, however, that some of the O1 sequences have unknown framework classifications. Altogether, we found a total of 16 cysteine frameworks associated with the designated major gene superfamilies (known and novel). Notably, all eight sequences classified under the B1 superfamily have no cysteine framework classifications.

A total of eight (one from C. eburneus and seven from C. tessulatus) predicted conopeptides classified as belonging to the B1 gene superfamily. Upon further analysis, the predicted mature peptide sequences from these transcripts were found to be highly similar with previously reported conantokin precursor peptide sequences ( Supplementary Table S3 ). Interestingly, one conantokin precursor in the database isolated from C. eburneus was not observed in the transcripts obtained from samples of the same species in this study, but was the top hit for three C. tessulatus B1 precursor peptides with over 90% precursor sequence similarity. Furthermore, alignment of the Tesseliconus B1 precursor sequences showed that although the signal and pro regions of the peptides are relatively conserved, the mature regions appear to be much more divergent ( Supplementary Figure S4 ).

Apart from the B1 superfamily conopeptides, we were also able to predict putative conopeptide precursors without known gene superfamily and cysteine framework classifications that have high amino acid sequence similarities with previously reported disulfide-poor conopeptides in the ConoServer database. In particular, we were able to identify a total of 54 putative disulfide-poor conopeptide precursors in C. eburneus and 56 in C. tessulatus. Table 4 shows a summary of the additional disulfide-poor conopeptide classes predicted from the two Tesseliconus species, while the individual sequence and ConoServer match details for all these putative conopeptides from C. eburneus and C. tessulatus are listed in Supplementary Tables S4 , S5 , respectively.

For both C. eburneus and C. tessulatus, most (>50%) of the mature peptide were found to have high sequence similarities with those obtained from species with known vermivorous diets, while the remaining peptides have top matches that are almost evenly split between those obtained from Conus species with piscivorous and molluscivorous diets ( Table 5 ). In terms of expression, conopeptides that are highly similar to those found in other vermivorous species registered the highest cumulative average expression in both Tesseliconus species ( Supplementary Table S6 ). However, we found that the cumulative average expression of piscivorous-like conopeptides is substantially higher in C. tessulatus than in C. eburneus, whereas the converse is true for molluscivorous-like conopeptides.

In both Tesseliconus species, the most highly expressed vermivorous-like conopeptides closely match conomarphin precursor sequences obtained from Conus betulinus. This is consistent with a previous proteomic study that identified a number of conomarphin peptides from C. eburneus exhibiting high degrees of similarity with those from C. betulinus (Itang et al., 2020). Some of the highly expressed molluscivorous-like conopeptides had significant sequence similarities with those obtained from Conus victoriae and Conus marmoreus, particularly peptides classified under the M gene superfamily with cysteine framework III. On the other hand, piscivorous-like conopeptides with relatively high average expression levels appeared to have significant similarities with those found in Conus magus (A superfamily), Conus catus (O1 and O2 superfamilies), Conus ermineus (cysteine framework XXII) and Conus californicus (divergent superfamily).

Agglomerative clustering of delta conopeptide sequences, believed to be one of the conopeptide families that enabled the shift from vermivorous to piscivorous behavior (Aman et al., 2015), showed that mature sequences from piscivores and vermivores indeed clustered more closely than their molluscivorous homologs ( Supplementary Figure S5A ). Interestingly, at the precursor peptide level, the molluscivorous and piscivorous clusters appear to be more closely related, with the vermivorous cluster forming an outgroup ( Supplementary Figure S5B ). However, it must be noted that there are only five representative molluscivorous delta sequences included in this analysis, three of which were observed from Conus textile.