The signal and mature regions of putative conotoxins were predicted using ConoPrec (Gao et al., 2017; Yao et al., 2019). After discarding those transcripts with disrupted cysteine frameworks in the mature regions, putative conotoxins with hits in BLASTX Results were translated according to the identified reading frames and were manually inspected. Most of them were assigned to known superfamilies based on their similarity to the highly conserved signal regions from the public ConoServer database. The total number of distinct precursors for C. striatus was 428, of which 405 were assigned to 11 known gene superfamilies (A, B1, M, O1, O2, O3, R, S, T, W, and Z), two cysteine-rich families (Conkunitzin and Con-ikot-ikot), and three hormone families (Conopressin/conophysin, Conoinsulin, and Prohormone-4, see more details in Table S3 ; Figure 1 ). Another set of 15 conotoxin precursor sequences could not be assigned to any known conotoxin superfamily (termed “unassigned”), and through alignment analysis, we assigned eight of them to the new superfamily (NFS-1 and NFS-2; Table S3 ).

A total of 14 known cysteine frameworks were also described. The maximum proportion of cysteine framework observed was VI/VII (C-C-CC-C-C), which accounted for 46.09% ( Figure 1 ). We also identified 106 proteins belonging to four different protein families ( Table S4 ), namely protein disulfide isomerase (PDI), conoporin, ferritin, and transcription factors, with potential roles in the processing and application of conotoxins. For example, 33 PDIs may be involved in the formation of disulfide bonds in conotoxins (Wang et al., 2017; Zamora-Bustillos et al., 2021), ten conoporins possibly enhance the hemolytic activity of venom (Hu et al., 2011; Terrat et al., 2012; Barghi et al., 2015), and one venom-related protein elevenin and 32 transcription factors are highly conserved and ubiquitous in other venomous animals (Safavi-Hemami et al., 2014; Jouiaei et al., 2015; Pardos-Blas et al., 2019).

The most diverse conotoxin precursor superfamilies, A, O1, and W, accounted for 77.56% of the total putative conotoxins. Among them, 149 (34.81%) belong to the O1 superfamily, 133 (31.07%) to the peculiar cysteine-poor W superfamily, and 50 (11.68%) to the A superfamily ( Figure 2 ). Interestingly, these data are highly consistent with the dominant position of the previously reported O1 superfamilies and VI/VII cysteine patterns (Pi et al., 2006; Walker et al., 2009; Himaya et al., 2021). Some of the inferred mature regions in the conotoxin precursors didn’t have any cysteine framework, such as the W and Z superfamilies that belonged to superfamilies exclusively formed by members without cysteines. We combined the conotoxin datasets with the total transcript dataset of C. striatus and named the predicted 428 conotoxins as Ts-001 to Ts-428 ( Table S3 ). Through the public ConoServer database, we collected these conotoxins in C. striatus previously identified by traditional methods and NGS technology for classifying into 19 known conotoxin precursor superfamilies, namely A, B1, B2, G2, H, I1, I2, J, M, N, O1, O2, O3, P, S, SF-mi2, T, W, and Z ( Figure 1 ) (Walker et al., 2009; Jagonia et al., 2019; Himaya et al., 2021). Compared with the previously reported conotoxin sequence of C. striatus, the amount of conotoxins found in this study is the sum of the previous amounts. These three superfamilies (R, Conoinsulin, and Prohormone-4) and three cysteine patterns (XVI, XXII, and XXIV) have been reported for the first time in the C. striatus. In addition, 10 of the 428 identified conotoxins have the same mature region sequences as previously reported in C. ermineus, C. striatus, and C. geographus ( Table S3 ), and they are classified into the A, O1, O2, T, and M superfamilies, Con-ikot-ikot, and Conkunitzin. With the technical reform and skilled application of cone snail transcriptome sequencing, more potential conotoxin sequences were excavated to help us establish a systematic conotoxin library.

After removal of those transcripts with duplication or truncated mature region sequences, we identified 272, 346, 203, and 247 conotoxin precursor sequences from the four transcriptome datasets of RS, VD, VG, and SG, respectively ( Tables S5 - S8 ). The comparative distribution of putative conotoxin sequences is summarized ( Figure 2 ). Interestingly, each organ had many unique conotoxin precursors, and 50.7% (187) of all inferred precursor peptide sequences were common to all transcriptomes, and 63.8% (235) of transcripts were shared between at least three specimens.

The 187 conotoxin sequences shared by four organs were assigned to 14 superfamilies: A, M, O1, O2, O3, R, W, Z, Con-ikot-ikot, Conkunitzin, New superfamily, Conopressin/conophysin, Conoinsulin, and Prohormone-4 ( Figure 2 ). The O1 and W superfamilies demonstrated abundant diversity in the four organs, while the Z and W superfamilies were evenly distributed ( Figure 2 ). The diversity of members of the different conotoxin precursor superfamilies was uniform across the four organs. In addition, 346 conotoxin precursors were identified from the venom duct transcriptome, accounting for 80.8% of all identified conotoxin sequences, indicating that the venom duct was the most important factory for the organizational production of conotoxins ( Table S6 ).

To investigate the transcription levels of conotoxins in each organ, clean reads were mapped back onto the de novo assembled unigenes, and the values of fragments per kilobase of transcript per million (FPKM) mapped fragments were calculated (Pardos-Blas et al., 2019). The top 20 conotoxins (with the highest FPKM values) were selected from each dataset for comparison to find the expression levels of the most abundant transcripts. We observed that transcription levels in the VD (FPKM values above 10,426.99) were generally higher than those in the RS, VG, and SG ( Table S9 ).

Figure 3 illustrates the differences in the spatial distribution of different conotoxin superfamilies among the RS, VD, VB, and SG datasets of C. striatus. The M, R, Z, W, and new superfamilies were not found in VD ( Figure 3 ). Although there were significant differences in the relative transcription levels of the conotoxin superfamilies, the A and O superfamilies had the highest transcription level in the four organs. Most precursors had widely varying FPKM between specimens, for example, TS-022 was common in the VD transcriptomes (FPKM value of 31,827.03, which decreased to 5,881.78, 154.69, and 2,036.95 in the RS, VG, and SG, respectively, Table S9 ). The TS-022 was highly similar (only two aa substitutions) to the κ-conotoxins SIVA from C. striatus, suggesting its κ-conotoxin character and potential voltage-gated ion channel activity (Wang et al., 2003; Santos et al., 2004). Ts-018, with a mature region sequence matching the previously identified alpha-conotoxin SI of the A superfamily, had the highest transcription level in the VD. Similarly, Ts-165 was the highest expressing Ca⁺ channel targeted omega-conotoxin SVIA of superfamily O1 in VD (FPKM value 31,827.03), but it was the lowest in RS, VG, and SG (FPKM values less than 1,591.44). Ts-147 was a highly expressed (Abalde et al., 2018) conotoxin identified in both RS and VD of C. striatus and was assigned to the Contryphans superfamily, and through other studies, it was also found to be dominant in C. ermineus (Abalde et al., 2018; Himaya et al., 2021). The Z and W superfamilies had high sequence diversity despite their relatively low expression levels. These data depicted that most of the identified conotoxins were synthesized at high levels in the venom duct, and some conotoxins were expressed in the venom gland and the salivary gland at low levels. Surprisingly, we also observed that the quantity and expression level of conotoxins in RS were the second highest relative to those in VD ( Figure 3 ).

Using both NGS and TGS to sequence C. striatus specimens, the transcriptome databases were generated after ConoSorter searching and manual editing. To annotate conotoxin coding sequences among the predicted unigenes, we constructed a local reference database of known conotoxins from the public ConoServer database and manually checked them using the ConoPrec (Gao et al., 2017; Yao et al., 2019). By contrast, the NGS and TGS were quite distinct. Only a few conotoxins were easily detectable with the NGS, while the TGS detected a more complex dataset. The TGS yielded 366 conotoxin precursors that were clustered into 11 known gene superfamilies (A, B1, M, O1, O2, O3, R, S, T, W, and Z) and two cysteine-rich families (Conkunitzin and Con-ikot-ikot; Tables 2 and S10 ), while the NGS identified fewer conotoxin precursors (82) that were classified into nine known gene superfamilies (A, B1, M, O1, O2, O3, R, T, and W) and two cysteine-rich families (Conkunitzin and Con-ikot-ikot; Tables 2 and S11 ). Most of these conotoxin sequences were reported for the first time, and some possess new cysteine patterns. These predicted conotoxins were summarized and named after the two technologies as Ts-NGS-001 to Ts-NGS-082 and Ts-TGS-001 to Ts-TGS-366 ( Tables S10 and S11 ).

These sequences from the same sample have only 4.7% overlapping between the TGS and NGS platforms, namely the A, B1, O, R, T, and new superfamilies, Conkunitzin and Con-ikot-ikot families, and three hormone families (Conopressin/conophysin, Cononsulin, Prohormone-4; Table S12 and Figure 4 ). Except for these overlapping conotoxins, 62 and 346 conotoxins were obtained from NGS and TGS, respectively. Among them, the A, O1, and W superfamilies have significant quantitative differences ( Figure 4 ). A total of 151 O1 superfamily conotoxins were identified, and among them, 141 were found by the TGS. Apart from the typical VI/VII cysteine framework, we also discovered I, XIV, XVI, and XVIII (NGS only VI/VII). A total of 133 W superfamily conotoxins were identified without non-cysteine frameworks, of which the TGS-derived conotoxins accounted for 79% ( Table 3 ). Meanwhile, 47 conotoxins belonging to superfamily A were found by TGS. Although their number was much larger than that of NGS, the cysteine patterns were relatively consistent. The S and cysteine-poor Z superfamilies were only discovered by TGS, even though their expression was low. The TGS found many new toxins and cysteine patterns that NGS did not identify, but their overall transcription levels were very low, especially for the A and W superfamilies ( Figure 4 ). The TGS dominates the O superfamilies with a much higher number and transcription of conotoxins than the NGS. The emergence of the TGS technology seems to be reliable in generating a more comprehensive library of conotoxins.

A total of 428 unique conotoxin precursors were classified into 11 known superfamilies, two new superfamilies (NSF-1 and NSF-2), and unassigned superfamilies based on the sequence similarity in C. striatus ( Table S3 ). A total of 216 typical conotoxin sequences from A, O, M, and W superfamilies and new superfamily were selected to construct a phylogenetic tree ( Figure 5 ). The phylogenetic evaluation displayed that the branches of these conotoxin superfamilies from the same Conus species appeared in different lineages. NSF-1 and NSF-2 superfamilies could not be assigned to any known families but were independently distributed among O1, O2, and O3 superfamilies ( Figure 5 ).

The A superfamily, as a typical conotoxin group in piscivorous Conus species with very distinct structural and functional diversity, can be subdivided into two subgroups: (i) alpha (α) with cysteine frameworks I and II (Azam and McIntosh, 2009), such as those that preferentially target nicotinic acetylcholine receptors (nAChRs) and ultimately inhibit neuromuscular transmission producing paralysis; and (ii) kappa (κ) with cysteine framework IV, which selectively targets voltage-gated potassium channels (K⁺) (Puillandre et al., 2010; Robinson and Norton, 2014). For the A superfamily, 50 conotoxins were determined in this study, and ten sequences have the common I pattern (CC-C-C framework), 11 sequences share the common II pattern (CCC-C-C-C), and 30 transcripts own the IV cysteine pattern (CC-C-C-C-C). This is consistent with the previously reported abundance and variety of A superfamily conotoxins in other piscivorous Conus species (Hu et al., 2011; Terrat et al., 2012; Safavi-Hemami et al., 2014). However, the typical paralytic Conus peptides found in fish-hunting cone snails are from α-conotoxin families. Among them, the representative α3/5 subfamily has the cysteine framework I. These peptides, such as Ts-018, usually have 3 amino acids between cys2 and cys3 and 5 amino acids between cys3 and cys4 ( Figure 6 ). Ts-014 targeting α1β1γδ nAChRs has the following IV cysteine pattern (Himaya et al., 2021). Ts-023 and Ts-032 have cysteine framework IV, which leads to characteristic spastic paralysis symptoms by inhibiting voltage-gated potassium channels (Craig et al., 1998). In addition, the arrangement of cysteine skeleton residues with the primary structure cysteine framework IV (i.e., CC-C-C-C-C-) belongs to the characteristic signs.

The M superfamily conotoxins are dominant in terms of transcription abundance and diversity of cysteine frameworks (XXXII, I, II, III, IV, VI/VII, IX, XIV, XVI) (Rajesh, 2015, Jacob and McDougal, 2010, Zhou et al., 2013). According to their pharmacological targets, such as voltage-gated ion channels (Na⁺, K⁺, and Ca²⁺) and ligand-gated ion channels (nAChR, AMPAR, GABAR, and 5-HT3R), these conotoxins are further grouped into different pharmacological families (α, ι, κ, μ) (Jacob and McDougal, 2010; Rajesh, 2015; Franco et al., 2018). A total of nine M superfamily conotoxins presented with distinct cysteine frameworks ( Figure 6 ), including four with III (CC-C-C-CC), two with II (CC-C-C-C), and one with XVI (C-C-CC) cysteine frameworks. μ-Conotoxins Ts-052 and Ts-058 have cysteine framework III, which belongs to S3-G04 and SIIIA conotoxin, respectively, act on Na⁺ channels, and have analgesic activity (Bulaj et al., 2005; Schroeder et al., 2008; Yao et al., 2008). In addition, the cysteine-free components of Ts-061 have been previously classified as Conomarphins, with only five sequences described in mollusk-hunting Conus species (Akondi et al., 2014).

The O superfamily conotoxins have high transcription levels, and several isoforms have been proved to play important roles in capturing prey and/or defense. In general, the O superfamily has been subdivided into O1, O2, and O3 and is further classified into δ, µ, ω, κ, and γ families with an extensive target repertoire in ion channels (Gilly et al., 2011; Dutertre et al., 2013; Gao et al., 2017). They are the most abundant group in terms of conotoxin number ( Table 3 ), with the extensive I, VI/VII, XIV, XVI, XVIII, and XXIV patterns of cysteine frameworks. Most O1, O2, and O3 superfamily members are unusually hydro-phobic peptides with the same cysteine framework of C-C-CC-C-C and with a similar inhibitory cystine knot motif (Jiang et al., 2017; Bernaldez-Sarabia et al., 2019). A total of 149 O1 superfamily sequences were identified in our present study. Among them, 103 presented the framework VI/VII. This typical isoform was a dominant family in quantity and displayed a high transcription level (such as an FPKM of 37,365.27 for TS-165) in the venom duct. TS-165 is ω-conotoxin SVIA that acts on presynaptic membranes ( Figure 6 ). It binds and blocks Ca²⁺ channels with the cysteine pattern of VI/VII (Nielsen et al., 1996). Fourteen O2 superfamily conotoxin precursors were identified, which were designated Ts-211 to Ts-224. Two unique O2 superfamily (cysteine framework VI/VII) and 12 Contryphan conotoxins were identified, while only one (Ts-213) had been identified previously ( Figure 6 ). Contryphans are short single disulfide-containing conotoxins (only 8 aa) that display diverse functions, generally targeting Ca²⁺ channels (Hansson et al., 2004; Sabareesh et al., 2006; Figueroa-Montiel et al., 2018). One O3 superfamily precursor was identified in C. striatus, with the classic signal peptide sequence ( Figure 6 ) and the cysteine-rich framework of VI/VII (Bernaldez-Sarabia et al., 2019; Dutt et al., 2019).

Conoinsulins regulate carbohydrate and fat metabolism and act as neuromodulators of energy homeostasis and cognition. Insulin is usually synthesized as a precursor with three regions (A, B, and C chains), from which proteolytic cleavage of the C peptide in the Golgi releases the insulin heterodimer with A and B chains connected by disulfide bonds (Ward and Lawrence, 2011). Some Conus species weaponize their endogenous Conoinsulins to induce insulin shock (decreases blood glucose levels) to capture prey. Three fish-hunting cone snails, such as C. ermineus (Abalde et al., 2018), C. geographus (Safavi-Hemami et al., 2014), and C. tulipa (Dutt et al., 2019), have evolved specialized insulins, which have much greater similarity to fish insulins than to the molluscan hormones (Safavi-Hemami et al., 2015). Here, we identified two putative venom Conoinsulins in C. striatus, which were highly similar to Con-ins Im1 and Con-ins Im2 from C. imperialis. (Himaya et al., 2021). Multiple sequence alignment of insulin sequences from humans, zebrafish, and cone snails depicted that the primary sequence of insulin family was poorly conserved except for the strictly conserved cysteines (Lu et al., 2020) ( Figure 7 ). By contrast, Conus venom Conoinsulins are greatly variable and differ from vertebrate insulins by having one additional disulfide bond between A and B chains. In brief, venom insulin sequences exhibit pronounced divergence, and the few relatively well-conserved amino acids include Gly5, IIe6, Glu19, and Tyr25 in the A chain, Thr4, Leu7, Pro16, Gly21, and Leu13 in the B chain, and all eight cysteines.

The W superfamily conotoxins, with non-cysteine mature peptides and original report in C. marmoreus, were particularly abundant in C. striatus and C. ermineus (Hansson et al., 2004; Dutertre et al., 2013; Pardos-Blas et al., 2019). A total of 133 new W superfamily sequences were identified in our present study, accounting for 31% of conotoxin precursors ( Figure 7 ). This W superfamily lacks cysteine mode, and its biological activity has not been systematically studied. Therefore, further analysis is needed to determine the pharmacological effects of the W superfamily for exploring the evolution of various conotoxin superfamilies.

Conosorter compared the conserved amino acid sequences in conotoxin sequences and found that the eight conotoxin precursors had sufficient differences from the sequences of any known superfamily. Therefore, the eight conotoxins were divided into NSF-1 and NSF-2 signal superfamilies through software analysis and phylogenetic analysis ( Figure 7C ). In the two new superfamily conotoxin sequences, there were only a signal region and a mature region, and the precursor region was lacking. We unexpectedly found that the number of cysteines in these sequences was singular, which may form dimers. The NSF-1 superfamily sequence contained a total of 90 aa and had a C-C-C cysteine pattern. The NSF-2 conotoxin precursors had a signal region sequence (MKAVAVFMIVALAVAYG) and an XXII disulfide framework, while the physiological and pharmacological properties of their mature peptides should be characterized with more functional studies (Pardos-Blas et al., 2019).

In addition to the above-mentioned major superfamilies, several less representative B1, R, S, T, and Z superfamilies and the Con-ikot-ikot and Conkunitzin families were discovered in the four transcriptomes. This study identified R, S, T, and Z superfamiles sequences with the cysteine-free pattern. Nine Conkunitzin family sequences presented with IX (C-C-C-C-C-C-), XXII (C-C-C-C-C-C-C-C-C), or VIII (C-C-C-C-C-C-C-C-C-C) pattern ( Table S3 ). Meanwhile, 12 Con-ikot-ikot peptides with the novel (C-CC-C-C-C-CC-C-C-C) framework were identified, which were reported to specifically target post-synaptic AMPA receptors (Walker et al., 2009). Interestingly, no transcriptomic evidence of I superfamily conotoxins was found, although they are often widespread among other piscivorous cone snails.

Previous studies have revealed that each complete conotoxin precursor has a hypervariable mature region, a fairly conserved precursor region, and a hyperconserved signal region (Wu et al., 2013; Yang and Zhou, 2020). However, we found that the diversity relationship among the mature, precursor and signal regions is more complex than reported. In the present study, some conotoxin precursors had the same hypervariable mature regions, but the putative hyperconserved signal regions were different due to point mutation [including insertions and deletions (indels) and/or substitutions].

All 172 identical mature regions of different conotoxins were aligned to detect point mutations and indels in the signal regions, among the A, M, O, and W superfamilies, these variances account for 25%, 1%, 58%, and 16%, respectively ( Table S13 , Figure S2 ). As listed in Table S13 , the conotoxin sequences have the same mature region; 74 and 71 conotoxins were deleted and replaced in the signal region, and 11 and 31 conotoxins were deleted and replaced in the precursor region, respectively. A total of 146 point mutations were present in the signal regions, suggesting that point mutations are frequent in the signal regions but infrequent in the precursor regions ( Figure 8 ). In addition, 158-point mutation sequences were identified by the TGS technology, and only 14 sequences were found by the NGS technology, indicating that the TGS broke the traditional view of the conservation of conotoxin signal peptides and the variability of mature peptides obtained by the traditional NGS technology.

Three C. striatus specimens, about 9 cm in length, were collected from the offshore areas of the West Island near Sanya, South China Sea, China. These cone snails were shipped to the indoor laboratory at Hainan Medical University (Haikou, Hainan, China), and then were kept in an ecological breeding system for about 3 days. Four organs of the cone snails, including radular sheath, venom duct, venom gland, and salivary gland, were dissected and then immediately stored at –80°C before use.

Total RNA was extracted from these samples using QIAzol reagent (Qiagen, Fredrick, MD, USA). A NanoDrop nucleic acid quantifier was used to detect the purity of RNA (A260: A230 ratio of 2.0–2.2; A260:A280 ratio of 1.8–2.1), a Qubit RNA assay (Invitrogen, Carlsbad, CA, USA) was used to accurately quantify the RNA concentrations, and an Agilent 2100 Bioanalyzer (Palo Alto, CA, USA) was used to detect the integrity of RNAs.

The high-quality RNA samples (OD260/280 = 1.8–2.2, OD260/230 ≥ 2.0, RIN ≥ 8, > 1 μg) were used to construct individual sequencing library. PolyA mRNA was purified from total RNA using oligo-dT-attached magnetic beads and then fragmented by a fragmentation buffer. Taking these short fragments as templates, first-stranded cDNA was synthesized using reverse transcriptase and random primers, followed by second-stranded cDNA synthesis. Then, the synthesized cDNA was subjected to end-repair, phosphorylation, and “A” base addition according to the library construction protocol. Afterward, sequencing adapters were added to both sides of the cDNA fragments. After PCR amplification of cDNA fragments, 150–250 bp targets were obtained. Finally, we performed the paired-end sequencing on an Illumina HiSeq X Ten platform (Illumina Inc., San Diego, CA, USA).

Total RNA was extracted by grinding tissues in TRIzol reagent (Cat# DP424, TIANGEN Biotech Co. Ltd, Beijing, China) on dry ice and was processed following the manufacturer’s instructions. The integrity of the RNA was determined with the Agilent 2100 Bioanalyzer and agarose gel (Lonza) electrophoresis. The purity and concentration of the RNA were determined with the Nanodrop micro-spectrophotometer (Thermo Fisher Scientific, USA). The mRNA was enriched by Oligo (dT) magnetic beads. Then, the enriched mRNA was reverse transcribed into cDNA using a Clontech SMARTer PCR cDNA Synthesis Kit (TAKARA). PCR cycle optimization was used to determine the optimal amplification cycle number for the downstream large-scale PCR reactions (PrimeSTAR GXL DNA polymerase). Then, the optimized cycle number was used to generate double-stranded cDNA. The large-scale PCR was performed for the next SMRTbell library construction. cDNAs were damage-repaired and end-repaired, and the sequencing adapters were ligated using a SMRTbell Template Prep Kit (PacBio). Each SMRT bell template was annealed to a sequencing primer, bound to polymerase, and sequenced on a PacBio sequel platform by Nextomics Biosciences Co. Ltd. (Wuhan, Hubei, China). Circular consensus sequences were obtained using the SMRTlink 5.1 software (PacBio). A total of 35.6 Gb subread sequences were obtained by correction between subreads. According to the existence of 5’-end primer, 3’-end primer, or polyA tail, these sequences were divided into FLNC reads and non-full-length sequences. The former was clustered by iterative clustering in Iterative Clustering for Error Correction (Wang et al., 2021) algorithm software to generate the cluster consensus sequences. Subsequently, we corrected the polished consensus sequences of the TGS data through LoRDEC software (PacBio) with default parameters, and any redundancy in corrected consensus reads was removed by CD-hit (version 4.7) (Li and Godzik, 2006) to obtain the final transcripts for subsequent analysis.

Illumina NGS data from the same samples were assembled with Trinity to produce unigenes (Shimizu et al., 2006). Transcripts were annotated by searching against five public databases, namely, clusters of orthologous/prokaryotic groups of proteins (Tatusov et al., 2001), NCBI NR (Li et al., 2002), GO (Ye et al., 2006), KEGG Ortholog (Kanehisa et al., 2017), and SwissProt (Soudy et al., 2020; UniProt, 2021). Transcripts >500 bp were selected for SSR analysis using the MIcroSAtellite identification tool (MISA; http://pgrc.ipk-gatersleben.de/misa, http://pgrc.ipk-gatersleben.de/misa/) (Thiel et al., 2003). MISA can identify seven SSR types, mononucleotide, dinucleotide, trinucleotide, tetranucleotide, pentanucleotide, hexanucleotide, and compound SSR, by analyzing transcript sequences. We used ANGEL software (Shimizu et al., 2006) (https://github.com/PacificBiosciences/ANGEL ) to perform CDS and protein predictions for all isoforms. Using a custom Python script (Wang et al., 2018) (alternative_splice.py, https://github.com/Nextomics/pipeline-for-isoseq ), AS events were identified from alignments and classified into exon skipping, intron retention, alternative donor site, alternative acceptor site, and alternative position. We used the ncRNA pipeline process (Li et al., 2014) (https://bitbucket.org/arrigonialberto/lncrnas-pipeline ) to obtain the LncRNA redundant transcript sequence.

Homologous searches and ab initio prediction methods (Kaas et al., 2012) were used to predict conotoxins from transcriptome datasets. All previously known conotoxins were downloaded from the ConoServer database for homologous prediction to build a local reference dataset. We subsequently used BLASTX (with an E-value of 1e-5) to run our assembled sequences against the local dataset. Those unigenes with the best hits in the BLASTX data were translated into conotoxins sequences. In addition, an ab initio prediction method using Hidden Markov Model (HMM) (Laht et al., 2012)was adopted to discover new conotoxins. The four conotoxin datasets were grouped into different classes according to the superfamily and family classification in the ConoServer database. Sequences of each familes were aligned with ClustalW, and the ambiguous results were checked by manual correction. Finally, a profile HMM was built for the conserved domain of each class using hmmbuild from the HMMER 3.0 package (http://hmmer.janelia.org ), and the hmmsearch tool was then applied to scan every assembled unigene/EST for identification of conotoxins (Wheeler and Eddy, 2013).

Conserver is a machine program for large-scale recognition and classification of conotoxins (Kaas et al., 2012). The Conoprec tool implemented in the Conserver database was used to identify the three domains (signal, propeptide, and mature) of the predicted conotoxin precursors and the cysteine frameworks of the mature functional peptides. Precursors were assigned to different protein superfamilies based on the percentage of sequence identity (> 70%) to the highly conserved signal region. If the conservation of the signal region was lower than the threshold of any reported conotoxin superfamily, members of the new superfamily were named “NSF-1 and NSF-2” according to the conservation of the conotoxin precursor peptide (Robinson and Norton, 2014). Those conotoxins without signal regions or incomplete signal regions but still showing similarity in either the pro- or mature region were classified as the “unassigned” group.

Multiple sequence alignments for all conotoxins were performed using MUSCLE (v. 3.8.31) (Edgar, 2004), followed by manual adjustment using the GeneDoc software. A phylogenetic analysis was carried out using the maximum likelihood method with RAxML8.1, and statistical supports were assessed using 1,000 bootstrap pseudo-replicates. We used the complete conotoxin to perform a phylogenetic analysis of the highly diverse A, M, O, W, and new conotoxin precursor superfamilies. The amino acid sequences were aligned using ClustalW 2.1 software (Thompson et al., 1994). We chose maximum parsimony as a tree-building method in MEGA 7.0.14 software (Kumar et al., 2016).

Clean reads were mapped back onto the assembled precursors, and transcript abundance was calculated in FPKM. The FPKM value of 0.1 was used as a threshold for judging whether the gene was expressed. RSEM software (v1.1.12) was used to calculate the expression levels, and the unigene set clustered by TGICL software was used to build a library (Pertea et al., 2003; Li and Dewey, 2011). Bowtie2 software was used to align the transcriptome sequence of each sample to the unigene library, and the transcription levels of all unigenes in different samples were counted (Langdon, 2015).

The presented multi-omics data were visualized using Prism software v7.0.0, and Venn maps were generated using Venny v2.1.0 software.