Sequencing as well as assembly of the transcriptome of Exaiptasia diaphana was developed and deposited at the National Center for Biotechnology Information (NCBI) database (BioProject: PRJNA694232, and SRA accession: SRP303227). High-throughput sequencing revealed about 22.31, 23.38 and 20.75 million paired-end short reads with a length of 150 bp and Q30 were 92.47%, 92.47% and 92.77% for small, middle and big sea anemones respectively ( Table 1 ). After filtration of low-quality reads as well as adapters, 21.94, 23.02 and 20.42 gigabases (Gb) of clean data were obtained for every sample. Clean data ratio was projected to be 98.34%, 98.46% and 98.41% respectively ( Table 1 ). For the whole transcriptome, De novo assembly of short reads produced 180210 unigenes, and N50 was determined to be 1271 bp with mean lengths of 794 bp. Average respective GC levels of contigs and unigenes for the transcriptome were 40.58% and 42.65% ( Table 1 ). To assess the functions of unigenes, annotations were done based on 4 databases. The 180210 unigenes were allocated into the databases, Uniprot (48716 unigenes), Nr (48631 unigenes), KOG (25973 unigenes), and KEGG (21364 unigenes), while 124166 more unigenes were not annotated in these databases ( Supplementary Figure 1 ).

About 27769 genes were enriched in 34 KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways, and were allocated to 5 primary categories: environmental information processing (3844), cellular processes (3892), genetic information processing (3804), organismal systems (6936) and metabolism (9293) ( Supplementary Figure 2 ). Most of the unigenes were allocated to metabolism, with the global as well as overview maps exhibiting the highest counts of annotated unigenes (3414). Another GO analysis illustrated that 139662 unigenes were annotated into 3 categories: cellular components (37080), biological processes (58038), and molecular functions (44504) with binding to be the most enriched one (18132), respectively ( Supplementary Figure 3 ). We also annotated the functions of the predicted proteins using Cluster of KOG database. A total of 23202 unigenes were allotted KOG annotations and were classified into 25 molecular families ( Supplementary Figure 4 ).

In this study, 533 identified transcripts were assigned into 75 groups based on predicted functions, which were projected based on sequence homologies at amino acid levels with peptides as well as proteins from UniprotKB databases, among which 72.98% were proteins and 27.02% were peptides ( Figure 1 and Supplementary Table 1 ). Majority of protein constituents corresponded to metalloproteases, chymotrypsinogen like, pancreatic lipase-related protein like (PLRP like), G-protein coupled receptor, and collagen, while peptide components corresponded to the ShK domain, thrombin, Kunitz-type, and insulin-like peptide and defensin ( Figure 1 and Supplementary Tables 2, 3 ). In general, a similar finding has been reported in transcriptomes of other sea anemones (Madio et al., 2017; Ramirez-Carreto et al., 2019).

A total of 166, 181, 148 putative proteins and peptide toxin precursors were respectively identified from Small, Middle, and Big datasets from three body-sized Exaiptasia diaphana. Comparative distributions of protein as well as peptide toxins is shown in Figure 2A . Of these putative proteins as well as toxin precursors, 31 were detected in all three datasets (Ed-S, Ed-M and Ed-B), 15 were common to both Ed-S and Ed-M, 16 were shared by Ed-S and Ed-B, and 24 were common to both Ed-M and Ed-B. The 31 common toxins precursors were classified into 18 superfamilies, such as CAP, ShK domain, metalloproteinase, triacylglycerol lipase, and defensin.

Fragments per kilobase of exon per million mapped fragment (FPKM) values were evaluated to denote transcription levels for every protein and peptide toxin. The top 10 protein as well as peptide toxins (with highest FPKM values) were chosen from every dataset ( Figure 2B ), and ShK domain, triacylglycerol lipase, CYP17, metalloproteinase, as well as C-type lectin were transcribed in 3 various developmental stages of Exaiptasia diaphana. The ShK domain dominantly expressed. Moreover, in the middle development stage of Exaiptasia diaphana, collagen levels were the highest, and were comparatively high in the mature stage.

Sea anemones are sessile animals, so peptide neurotoxins are very important in immobilization of prey as well as defenses against predators. Peptide neurotoxins are among the best characterized constituents of sea anemone venoms in terms of their action mechanisms (Madio et al., 2019). They associate with various ion channels, such as VGSCs (Bosmans and Tytgat, 2007), VGPCs (Castaneda and Harvey, 2009), ASICs (Cristofori-Armstrong and Rash, 2017), and TRP channels (Monastyrnaya et al., 2016; Logashina et al., 2017). Here, 144 peptide toxins were obtained from the Small, Middle, Big, and Combine datasets of Exaiptasia diaphana, and named as Ed-001 to 144 in order ( Supplementary Table 3 ). Cysteine-rich peptides are present in almost all sea anemone neurotoxins, and diversities of structural scaffolds among them is notable. The cysteine patterns of the 144 peptides were divided into 11 major categories, each with different subclasses ( Figure 3 ). The peptides containing 4 cysteines forming 2 disulfide bond patterns (IV) and 6 cysteines forming 3 disulfide bond patterns (VI) are the most abundant, and their numbers are 37 and 44, respectively. In addition to peptides with even number of cysteines, there are also a large number of peptides with odd number of cysteines.

So far, the 3D structure and/or cysteine-pattern of nine typical anemone peptide toxins have been identified, which are ShK, Kunitz-domain, β-defensin-like, Anemonia sulcata toxin III (ATX-III), boundless β-hairpin (BBH), inhibitor cystine-knot (ICK), epidermal growth factor-like (EGF-like), proline-hinged asymmetric β-hairpin (PHAB), and small cysteine-rich peptides (SCRiPs) respectively (Madio et al., 2019). The putative peptide toxins with typical and unique homologues were ShK domain (11 homologues), Kunitz-type peptides (10 homologues), β-defensin-like (5 homologues), and EGF-like (3 homologues).

The ShK domain was named after the ShK toxin, which had been identified from the venom of Stichodactyla helianthus (Castaneda et al., 1995). This toxin is a member of Kv type 1 family, which blocks Kv1.1, Kv1.2, Kv1.3, Kv1.6, Kv3.2, and Kca3.1 channels (Kalman et al., 1998; Rauer et al., 1999; Yan et al., 2005). ShK blocks Kv1.3, which has important functions in T and B lymphocyte subsets involved in autoimmune conditions. Therefore, ShK is potential immune modulator for treatment of autoimmune diseases (Pennington et al., 1995). ShK analogs have been developed to be Kv1.3-specific (Beeton et al., 2011). The ShK-186 analogue is being developed as a treatment for autoimmune diseases (Chi et al., 2012). BgK, an ShK domain peptide from Bunodosoma granulifera, a sea anemone, blocks Kv1.1, Kv1.2, as well as Kv1.3 potassium channels (Cotton et al., 1997). Here, eleven homologous sequences with ShK and BgK were identified, and their cysteine pattern was C-C-C-CX₃CX₂C, including the connection mode of disulfide bond as C1-C6, C2-C4, C3-C5 ( Figure 4A ). Among these sequences, by using blast alignment in NCBI, Ed-079, Ed-080 and Ed-081 are the same as previously reported sequences (GenBank No. KXJ15362.1, XP_020909625.1 and KXJ06239.1) from Exaiptasia diaphana in the database, respectively. The homology comparison results show that the sequence identity of Ed-079 and ShK is 43.75%, and the sequence identity of Ed-081 and BgK is 33.33%. The homology modeling prediction results show that Ed-079 and Ed-81 have similar 3D structures, indicating that they may also act on potassium channels ( Figure 4B ). As analogues of ShK and BgK, they are worthy of in-depth study in the future.

Epidermal growth factor (EGF) family-associated proteins that have been identified in invertebrates. Among them, L-EGF is a growth factor secreted from Lymnaea stagnalis, a gastropod mollusk. Gigantoxin I (ω-stichotoxin-Sgt1a), a peptide toxin from Stichodactyla giganteus, a sea anemone, can paralyze crabs (Hermann et al., 2004). With regards to sequence homology, this toxin has EGF activities as shown by rounding of human epidermoid carcinoma A431 cells (Shiomi et al., 2003). Moreover, gigantoxin I modulates the activities of TRPV1 channels (Cuypers et al., 2011), however, these activities result from involvement of the EGF receptor/PLA2/arachidonic acid/lipoxygenase pathway in indirect TRPV1 activation. Gigantoxin I was the first toxin that initiated this effect. Therefore, these toxins may elucidate on TRPV1 channel regulation (Cuypers et al., 2011). Three homologous sequences with Gigantoxin I were identified, and their cysteine pattern was C-C-C-CXC-C ( Figure 5A ). The homology comparison results show that the sequence identity of Gigantoxin I, Ed-94, Ed-95 and Human EGF are 35.71%, 36.36% and 34.09%, respectively. Taking Human EGF as the template, the homology modeling prediction results show that Gigantoxin I, Ed-93,Ed-94, and Ed-95 have similar three-dimensional structures, respectively ( Figure 5B ). Given that mammalian EGF as well as its family members are regularly associated with neoplastic diseases, the invertebrate EGF family members may inform the designs of erbB receptor antagonists.

Kunitz-type peptides have been discovered from venoms as well as glands of various venomous animals, where they are involved in protease inhibition and ion channel modulation (Mishra, 2020). In sea anemones, Kunitz-type peptides act on TRPV and type II Kv channels (Monastyrnaya et al., 2016; Gladkikh et al., 2020). Various sea anemone type II Kv channel toxins have dual activities, where they acts as blockers of potassium channels as well as serine protease inhibitors, which allows them to serve as defense molecules as well as neurotoxins for prey immobilization (Honma et al., 2008; Mishra, 2020). The AsKC1 and ShPI-I are some of the toxins that suppress proteases and Kv channels. The AsKC1 and ShPI-I belong to Kunitz-type peptide, which are derived from Anemonia sulcata and Stichodactyla haddoni, respectively (Schweitz et al., 1995; Garcia-Fernandez et al., 2016). Seven homologous sequences with ShPI-I and AsKC-1 were identified, and their cysteine pattern was CX₈CX₁₅CX₇CX₁₂CX₃C, including 3 disulfide bridges with C1-C6, C3-C5, C2-C4, connectivities ( Figure 6A ). Ed-100 was 98.28% similar to the thyroglobulin (GenBank XP_020905524.1) from Exaiptasia diaphana in the NCBI database. The homology comparison results show that the sequence identity of Ed-100, Ed-097 and SHPI-1 is both 51.85%. Taking ShPI-I as the template, the homology modeling prediction results show that AsKC-1, Ed-100, and Ed-097 are described by an α/β/α motif with a compact hydrophobic core ( Figure 6B ). To date, several Kunitz-type peptides have been retrieved from sea anemone’s, tentacles, secrete mucus, whole bodies and aggressive organs, including acrorhagi in various species of the Actiniidae family (Minagawa et al., 2008).

β-Defensins are essential antimicrobial peptides secreted as innate immune system constituents in various organisms (Suarez-Carmona et al., 2015; Shafee et al., 2016). But, in the sea anemone venom, β-defensin-like peptides are potential neurotoxins that can block ligand- and voltage-gated ion channels including Na_V types 1, 2 and 4, K_V type 3, and ASIC (Cariello et al., 1989; Ishida et al., 1997; Diochot et al., 1998; Diochot et al., 2004). APETx1 and CgNa belong to β-Defensin-like peptids. APETx1 is a 42-amino-acid peptide toxin of the sea anemone Anthopleura elegantissima, which suppresses voltage-dependent activation of hERG (Matsumura et al., 2021). CgNa was purified from Condylactis gigantea, sea anemone. It inhibits Nav types I and II (Salceda et al., 2007). Five homologous sequences with sea anemone toxin APETx1 and CgNa were identified, and their cysteine pattern was CXC-C-C-CC ( Figure 7A ). These cysteines are paired in a C1–C5, C2–C4, and C3–C6 fashion to form disulfide bonds that are importance for maintenance of the compact core β-defensins configuration. The homology comparison results show that the sequence identity of Ed-103 and APETx1, Ed-105 and CgNa are 37.84% and 36.59%, respectively. Taking APETx1 and CgNa as the template, the homology modeling prediction results show that Ed-103 and APETx1, Ed-105 and CgNa have similar 3D structures, respectively ( Figure 7B ).

The insulin-like peptides (ILPs) have been extensively characterized in invertebrates and vertebrates, with their genes constituting a bug multi-gene family, including genes that encode peptides in fishes, cone snails, amphibians, and mammals (Conlon et al., 1998; Caruso and Sheridan, 2011; Viscarra et al., 2013; Ahorukomeye et al., 2019). In vertebrates, the physiological role and structure of insulin are highly conserved (Conlon, 2001). In contrast, invertebrate insulins are highly variable with several functions, such as haemolymph glucose level regulation, neuronal signaling, reproduction, memory, and growth (Smit et al., 1998; Safavi-Hemami et al., 2016). In general, a majority of these ILP genes encode a single pre-propeptide with a signal peptide and adjacent B, C, and A peptides. In the process of peptide processing, signal peptides are first removed to form the main active proinsulin. Next, the proinsulin is transformed into mature peptide insulin by forming disulfide bonds between cysteine in A peptide and B peptide, and the C peptide was then removed by a proteolytic enzyme (Thompson and Di Gregorio, 2015).

Here, five ILPs (ILP-Ap01~05) were obtained from a transcriptome data of Exaiptasia diaphana, and the alignment sequence is shown in Figure 8 and Supplementary Table 1 . Multiple sequence alignment of insulin family members, involving cone snails, fishes, amphibians and mammals, showed that primary sequences of insulin family members were poorly conserved, particularly among invertebrates, apart from the stringently conserved cysteines (Lu et al., 2020). All ILPs from sea anemones exhibited classic cysteine characteristics shared by the entire protein family, while sequence identities for other amino acids with ILPs of other species were low. Two more cysteines were found in B chain-like domain as a phylum-specific character for ILPs in platyhelminths, which markedly differed from findings in other phyla. In addition, compared with human insulin, a short additional propeptide between B chain and signal peptides were predicted in all sea anemone ILPs ( Figure 8 ). ILP-Ap01, ILP-Ap02 and ILP-Ap03 are similar to human insulin and have contiguous B, C and A peptides. However, A-chain was inserted between pro-peptide and B-chain in ILP-Ap04 and ILP-Ap05, resulting in contiguous A, B and C peptides, even the IPL-Ap05 had no C peptide. Therefore, these novel findings indicate that the C chain is not necessary for the formation of insulin.

The homologous map of six ILPs shows a basic insulin structure, including residues that maintain an insulin-type hydrophobic core and consisting of A and B chains covalently linked by one intrachain disulfide bond and two interchain disulfide bonds (Mao et al., 2019). At the same time, the conserved cysteines with the characteristic amino acid sequence CCX₃CX₈C was also shown, which was the result of the retention distance in A peptide (Mayer et al., 2007). The projected tertiary structures for ILPs from sea anemone are shown in Figure 9 . The models established by the SWISS-MODLE (Waterhouse et al., 2018) and I-TASSER server (Zheng et al., 2019) for every ILP are presented only with the model exhibiting the highest C-core. The predicted tertiary ILP structures exhibited great similarities to those of proinsulin from human. The tertiary ILP structure was assigned into 3 main α-helices, and a C-domain loop separating helices 1 and 2. The typical hydrophobic core, which had been formed by the 3 α-helices was present in all ILP tertiary structures modeled by SWISS-MODEL and I-TASSER. However, all of the projected models exhibited same conserved helices as well as functional hydrophobic core-features that were present in all insulin family members.

Phylogenetic analyses of ILPs from different species are shown in Figure 10 . All ILPs from different species (sea anemone, fish, Conus, mammals, and amphibians) were clearly divided into five groups, and sea anemone ILPs had close evolutionary associations with Conus ILPs, relative to vertebrates ILPs.

Exaiptasia diaphana were collected in offshore areas of Lingshui City (18°40 E, 109°95 N), Hainan province, China, on August 17, 2020. A total of 3 Exaiptasia diaphana were collected at each developmental stage (including early stage, middle stage and mature stage). The mixture of 3 anemone whole body was used for total RNA extraction, respectively. RNA integrity number (RIN) values of extracted RNA were assessed by an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA). Then, 3 Illumina cDNA libraries were separately build using qualified RNA from various developmental stages. The libraries were sequenced on an Illumina HiSeq4000 platform (Illumina, San Diego, CA, USA) at BGI-Tech (BGI, Shenzhen, Guangdong, China).

SOAPnuke software (Li et al., 2008) was employed to filter reads with sequencing adaptors, as well as reads with > 10% of non-sequenced bases or > 50% of low-quality bases (≤10 was the base quality score). Trinity v2.5.1 (Haas et al., 2013) was used for assembling the remaining clean reads into contigs, which were then clustered based on sequence similarities and thereafter assembled into consensus unigenes via TGICL v2.1 (Pertea et al., 2003). We used the a Bowtie 2 read aligner to align clean reads with the de novo assembled assembly to evaluate values of gene transcriptions in assembled transcriptomes (Langmead and Salzberg, 2012). The established alignments were analyzed by RNA-Seq by Expectation Maximization (RSEM) v1.2.31 for transcript abundance estimation in FPKM terms (Li and Dewey, 2011).

This objective was achieved by searching various public databases (with E-value ≤ 10−5 as the threshold) including InterPro (Finn et al., 2017), Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa and Goto, 2000), NCBI Nr and Nt, UniProtKB/Swiss-Prot (Boeckmann et al., 2003), Clusters of Orthologous Groups (COG) (Tatusov et al., 2000) and Gene Ontology (GO). Blast2GO v4.1 (Conesa et al., 2005) was used to conduct GO annotations of Nr blast findings.

Homolog searches as well as an ab initio prediction approach was used for predicting sea anemone protein and peptide toxins from 3 transcriptome datasets. With regards to homologous predictions, sea anemone protein and peptide toxins from the BLAST database were used to build a local reference dataset. Then, BLASTX (with an E-value of 1e-5) was used to run the assembled sequences against a local dataset. Unigenes with best hits in BLASTX data were transformed into aa sequences. The three datasets of sea anemone protein and peptide toxins were allocated into various classes based on superfamily as well as family classification in BLAST database.

Manual inspection of predicted sea anemone protein and peptide transcripts was done by the BLAST database. Transcripts with duplicate or truncated mature region sequences were eliminated. Gene superfamilies, signal peptides, as well as cysteine frameworks of the peptide and protein toxins were checked for verification. Based on 75% identity of highly conserved signal peptide sequences, protein and peptide toxins were allocated into known superfamilies in the BLAST database. Those protein and peptide toxins with low similarity were classified as unknown groups.

Concatenated amino-acid alignments of the signal domains of all superfamilies were performed using MUSCLE v.3.8.3 (Edgar, 2004). Before phylogenetic analyses, divergent and indistinctly aligned blocks were eliminated by the Gblocks software. Then, JModelTest and Prottest3.2 (Darriba et al., 2011) were used for selection of the best-fitted model. Phylogenetic analyses were performed by the maximum likelihood approach with RAxML8.1, and the trees was visualized and annotated using the tree viewer of MEGA4 (Stamatakis, 2006). The 1,000 bootstrap pseudoreplicates were used to assess statistical supports.

The 3D structure model of putative peptides were developed from primary amino acid sequences via homology modelling on a SWISS-Model server in combination with threading-based LOMETS on an I-TASSER server (Waterhouse et al., 2018; Zheng et al., 2019). In the SWISS-Model server, the homologous sequences with sequence identify than 50% were selected as the template, and the QMEAN score was used to evaluate the model quality.

Sequence alignments as well as phylogenetic analyses of members of insulin/relaxin family from sea anemone, fish, Conus, mammals, and amphibians were acquired from UniProt ( Supplementary Table 3 ; www.uniprot.org/);. Boundary identifications for A and B chains of the peptides were conducted according to sequence annotations in this database. Sequence alignments were done by CLUSTALW, after which they were manually checked using the BIOEDI software. Phylogenetic analyses of platyhelminth ILPs were conducted using NJ and maximum likelihood approaches implemented in MEGA. The settings were: Poisson substitution models and bootstrap test approached with 1000 replications.