Author Correction: Proteomics Studies on the three Larval Stages of Development and Metamorphosis of Babylonia areolata

A correction to this article has been published and is linked from the HTML and PDF versions of this paper. The error has been fixed in the paper.


Results
Three Stages of Larval Development. Three groups of larval B. areolata at different developmental stages including the middle veliger stage before attachment (ZRZ-III), later veliger stage (velum atrophy) (ZRZ-V), and juvenile stage (ZRZ-VI) were used. Each group was replicated three times. The three stages of larval development are shown in Fig. 1a-c. The early veliger stage is distinguished by the appearance of the velum, feet, and shell. The middle veliger stage is a pelagic life for 1-2 d and the average shell height of the middle veliger larvae is about 520 μm, and the late veliger stage is characterized by a more rapid heartbeat and larger and elongated foot (Fig. 1a). After 10-12 d, the late veliger stage sinks to the substratum and leads a benthic life, the velum gradually degenerates and atrophies, after which the veliger completes its metamorphosis (Fig. 1b) and enters the juvenile snail stage (Fig. 1c).
Sodium dodecyl sulfate polyacrylamide gel electrophoresis of the extracted proteins. The determined protein concentration of the sample is shown in Table 1. The extracted proteins were analyzed by sodium dodecyl sulfate polyacrylamide gel electrophoresis (SDS-PAGE). The samples at the three developmental stages showed clear protein bands (Fig. 2), three parallel experiments were performed with good repeatability for each developmental stage.
Analysis of the Significant Differences in Protein Abundance. All of the raw liquid chromatographytandem mass spectrometry (LC-MS/MS) data from the three developmental stages were analyzed by one-way analysis of variance (ANOVA). Subsequently, the DEPs (p < 0.05) were selected. From the raw data, a total of 5,583 proteins, including 1,419 DEPs that satisfied the screening criteria, were present (ratio > ±2 and p < 0.05). A total of 721 protein groups were identified and 1,930 specific polypeptides were identified. The quantitative protein sequence information was extracted in batch from the self-built library. The resulting protein sequences were analyzed by BLAST, thereby obtaining a total of 1132 highly similar proteins. These proteins were ranked by similarity, and the top 28 proteins are shown in Table 2. The pairwise comparison of the three different developmental stages is presented in Table 3. The number of DEPs between ZRZ_VI and ZRZ_III was large, and 733 proteins were qualitatively different. Among the proteins with varying expression levels, 281 and 259 proteins GO Functional Annotations. GO describes genes in the organism and the attributes of gene products from three aspects, namely biological processes (BPs), molecular functions (MFs), and cellular component (CCs). The identified proteins were analyzed by GO annotations in three larval stages, the GO functional classification is shown in Fig. 4. Using the BP term, the number of proteins participating in metabolism and those involved in cellular and single organism processes was the first and second largest groups, respectively. Using the MF term, the number of proteins involved in binding and related to catalytic activity was the first and second largest groups, respectively. Using the CC term, the number of proteins involved in cells was the first largest, followed by those related to organelles, the macromolecular complex, and the membrane, respectively. Proteins with high similarity were annotated with GO terms, and the MF results are shown in Table 4. The results showed that most of these DEPs exhibited binding functions.

GO Enrichment Analysis of DEPs.
The GO annotation of the target protein set allowed the classification of these proteins according to BP, MF, and CC. The proportion of proteins in each class may indicate the effects of each of the three different developmental stages on each GO class. The differentially expressed proteins were analyzed by GO annotations, statistical analysis of the significantly enriched GO terms from data of the three developmental stages is shown in Fig. 5. The results showed that the most significantly enriched GO term was MF, organic substance metabolic process, and nitrogen compound metabolic process, in that order.

KEGG Pathway Enrichment Analysis of DEPs.
The KEGG pathway is used to analyze the significance level of each protein enriched in pathways, thereby determining proteins that significantly affect metabolism and signal transduction. The differentially expressed proteins were analyzed by KEGG in three larval stages, statistical assessment of the significantly enriched KEGG pathways is shown in Fig. 6. The results showed that the most significantly enriched KEGG pathways were the pathways related to ribosome and carbon metabolism, lysosome, focal adhesion, amino acid biosynthesis, the phosphatidylinositide 3-kinase/serine-threonine kinase (PI3K)/Akt signaling pathway, Epstein-Barr virus infection, purine metabolism, phagosome, glycolysis/gluconeogenesis, glyoxylate and dicarboxylate metabolism, and the tricarboxylic citrate cycle, in that order. The relationship between KEGG pathway and the number of protein sequence was showed in Fig. 7 and Table 5. The results showed that the ribosome pathway is the pathway that is annotated the most to protein sequence, carbon metabolism and lysosome successively.

Discussion
This study mainly focused on the three developmental stages of B. areolata, namely the veliger stage of larval B. areolata before and after attachment, and the juvenile stage. When entering the veliger prior to settlement stage (described as metamorphosing larvae), with the loss of velum and gill development, the larvae crawled with the developed foot and began their feeding. When morphological changes occurred, including the transition stage from planktonic larvae to benthic life and from the later veliger stage to the juvenile stage, the feeding habits of B. areolata also gradually changed from phytophagy to sarcophagy. Gastropods metamorphosis is important to understand the life and evolution of the species. In this study, we identified over a thousand proteins differentially expressed before and after metamorphosis. The results revealed two thirds genes of 17 differentially expressed proteins exhibited the same trends at the mRNA level certification

Proteins Related to the Metabolism in GO Functional Annotations.
In the BP-classified enrichment, the proportion of proteins related to metabolism was the largest (Fig. 4). According to GO functional annotations, the number of proteins involved in the metabolism of these molecules including dynein, ATPase, spectrin, and hemocyanin was the largest. These proteins activate important metabolic and cellular processes of larval B. areolata at the early development stage, thereby allowing the efficient and orderly development and metamorphosis of the larvae. Dynein is a high molecular weight protein that is distributed in nervous tissues and organs that function in the transport of metabolites 34 and may transfer vesicles along the microtubules, showing its important role in intracellular vesicular transport 35 . Dynein itself may not enable the transfer of vesicles or organelles, but may function with ATPase and dynactin 36 . Furthermore, dynein may regulate the proliferation and cycle of nerve cells, thereby promoting development of the nervous system 37 . There have been studies on distribution in the optic nerve system of molluscs 38 , the left and right spiral formations 39 , and reproductive development 28 . Spectrin, a cytoskeletal protein, may maintain the stability of membrane and its shape and participate in a variety of MFs 40 , and may also interact with actin to establish a membrane network that maintains the elasticity of the cell 41 . In addition, because spectrin is a multifunctional protein that plays a multifaceted role, it may participate in signal transduction, cell cycle, intracellular transport, immunological reaction, and other MFs. H. rufescens 42 , A. californica 43 and other molluscs have also been investigated. Hemocyanin, a respiratory protein, mainly transports oxygen to some invertebrates 44 . This protein is a complex structure in the hemolymph of gastropod and is essential for the circulatory and respiratory systems of gastropods. Hemocyanin also plays  Table 3. The number of differential proteins.  45,46 . Thus, this protein participates in phagocytosis, agglutination, and other immunological reactions 47 . H. diversicolor 48 and Rapana thomasiana 49 have also been studied.

KEGG Pathway Enrichment Analysis of Differentially Expressed Proteins in Three Stages.
In the KEGG pathway enrichment analysis, significantly enriched KEGG pathways were related to ribosome, carbon metabolism, and lysosomal pathways ( Fig. 6), indicating that at the developmental stage of larval B. areolata, the genes and proteins participating in these related pathways were considerably expressed and regulated.
Ribosome pathways. The ribosome is composed of rRNA and protein 50 , which are very important in the cell. Ribosomes are at the center of protein synthesis, and their levels of cells with quick proliferation and strong secretion is higher than that in the normal cell 51 . In eukaryotic cells, ribosome synthesis is a complex dynamic process in which hundreds of diverse genes are involved 52 . The results of this study showed that during the early development and metamorphosis of larval B. areolata, the genes and proteins that participate in the ribosome pathway are active and the ribosome is closely related to intracellular protein synthesis. This observation proved that in the early developmental stage, larval B. areolata may regulate various physiological and biochemical processes related to protein synthesis by regulating the genes participating in the ribosome pathway. Studies on Papana venosa 53 , H. diversicolor 54 , and some other gastropods have also been reported.
Carbon metabolism pathways. Carbon metabolism is the most basic of life activities, as it is an important pathway for energy metabolism in organisms. This pathway involves glycolysis, phosphopentose pathway, citric acid cycle, six carbon fixation processes, and methanol metabolism 55 . Carbohydrates provide energy and raw materials for embryonic development 56 . In this study, the signaling pathway related to glycolysis or gluconeogenesis was significantly enriched. In addition, many tissues and organs are formed and constructed during development; the fact that enrichment of the signaling pathway, in which DEPs participated, is related to the synthesis of amino acid also illustrated this point. The signaling pathway related to the digestion and absorption of protein was also significantly enriched, probably due to transition of the larval B. areolata from phytophagy to sarcophagy.
Lysosomal pathways. Lysosomes are membrane-coated cystic cells that contain a broad range of hydrolases, which can break down a variety of macromolecules and play significant roles in apoptosis and the cellular Pathways related to cytoskeleton, cell adhesion and ECM. "Skeletons" exist inside and outside of cells, including intracellular, transmembrane and extracellular components, "Skeletons" mainly refer to cytoskeleton, cell adhesion and ECM, respectively, and comprise an important network. KEGG analysis showed that focal adhesion was significantly enriched. Focal adhesion is a type of adhesive contact between the cell and extracellular matrix (ECM), which are mediated by integrins, resulting in dynamic cell anchoring connections 60 . The assembly and depolymerization of focal adhesions affect the cell shape and site and direction of pseudopods. In addition, extracellular and intracellular signal synergies regulate the formation, distribution, and turnover of focal adhesion spatially and temporally, which consequently affect cell migration. During cell migration, focal adhesion is a constant dynamic structure in the continuous links of assembly, depolymerization, and reassembly. Integrin, acting as a bridge, is a small and unstable adhesive complex that links ECM and intracellular protein formation 61 . When a large number of proteins are attracted to the adhesive complex, large and stable focal adhesions are formed. Cell mitosis is related to cell proliferation and differentiation, embryonic development, and tissue and organ formation. A series of processes such as integrin-mediated extracellular adhesion, focal adhesion assembly, stress fiber arrangement, and mitotic spindle orientation affect the mitosis of adherent mammalian cells. Myosin and dynein play important roles in these processes. In this study, the signaling pathway related to focal adhesion was significantly enriched and closely related to cell division in embryonic development. Moreover, myosin, dynein, and other proteins were identified. ECM receptors interact with the relevant pathway that is significantly enriched. The differentially expressed proteins are involved in functions of cytoskeleton, cell adhesion, and ECM in three    ECMs are composed of a complex mixture of proteoglycan, glycoprotein, aminoglycan, and other biomolecules; can interact with the aid of their surface receptors; and play important roles in cell adhesion, migration, differentiation, and proliferation 62 . The ECMs have potential roles in tissue remodeling, ECM remodeling has been proved essential in metamorphosis of amphibian 63,64 , insects 65 , and bivalve C. gigas 13 . This study result speculate that ECMs may function in regulating cell fates in metamorphosis of B. areolata. In this study, the pathway that related to the interaction of ECM receptors may complement each other with that related with focal adhesion and bridge cell adhesions well. The PI3K/Akt signaling pathway, which regulates cell metabolism and participates in the survival and antiapoptosis of cells, was also enriched. This pathway also plays an important role in the early embryonic development of mammalians 66 , and regulates embryonic stem cells 67 . Therefore, enrichment of this signaling pathway may considerably affect the development process of B. areolata.

Cluster Analysis of the Three Larval Stages of Development and Metamorphosis. Cluster anal-
ysis results showed that in the horizontal direction, the protein expression modes at ZRZ_III and ZRZ_V were of the same branch, and that the expression mode at ZRZ_VI was of another branch at the three stages of early development of larval B. areolata. This result is consistent with the developmental characteristics of larval B.   areolata and is more similar to those at ZRZ_III and ZRZ_V of the veliger stage than that at ZRZ_VI. Irreversible metamorphosis is necessary from the veliger stage to the juvenile stage; each larval undergoes a certain degree of change, which is largely reflected in the protein expression mode. Therefore, the cluster analysis results were consistent with the larval developmental process. In the vertical direction, DEPs were divided into two major branches. The protein expression in the first branch was upregulated with developmental stage. In addition, the proteins mainly participated in the formation of actin filament and corresponding biological process including cell adhesion, cytoskeleton construction, and ion transport. Actin is closely related with the growth of molluscs, and its genes are common in transcriptomics studies on the early development of molluscs. It mainly participates in actin filament formation, cell proliferation, intracellular signal transduction, as well as some other processes 68 .
Proteins that participate in cell adhesion play important roles in important physiological processes, including cell migration, growth, stress tolerance, and immunity 69 . Thus, the enrichment of each protein involved in physiological activities such as the cell migration and growth exhibited a down-and-up trend. The result showed that at the veliger stage, the physiological activities such as the cell migration and growth of B. areolata larvae were slightly weaker than those from the metamorphosis stage to the juvenile stage. The protein expression in the second branch was downregulated, and the proteins mainly participated in gene expression, protein synthesis and modification, and ribosome biosynthesis. In the development process from the early veliger to late stage, various parts of the organs of B. areolata larvae grow quickly. Consequently, the expression of proteins related to cell growth was upregulated. Ribosomal proteins were highly enriched in the B. areolata larvae in the early stage of development, as they participate in the synthesis and modification of proteins 70 that satisfy the needs of the rapid

Proteins Associated with Shell Formation, Body Torsion, Changes in the Feeding Habits, Attachment and Metamorphosis, and Immune Responses in B. areolata Larvae. Proteins related
to changes in feeding habits. In this study, a wide range of lipases, amylases, and proteases were identified. The digestive enzymes were diverse, consistent with the fact that B. areolata at different stages use different foods as a source of nutrition. The expression levels of these digestive enzymes varied at different developmental stages, showing the importance of regulating protein levels during adaption to changes in the feeding habits. In terms of the serine proteases, represented by trypsin and chymotrypsin, studies on crustaceans have shown that trypsin and chymotrypsin are the two most active proteases in the digestive gland, which account for more than 60% of the protein digestion in prawns 73 . In this study, the levels of chymotrypsin were in the the order of ZRZ_ VI > ZRZ_V > ZRZ_III, and same trend was observed at the mRNA level (Fig. 10). This is consistent with previous study that a chymotrypsin gene was highly induced during the metamorphosis of H. rufescens in the digestive system 74 and larva metamorphosis of C. gigas 13 , these results may show a transition of the digestive system in metamorphosis. Additionally, there were different levels of carboxypeptidase and serine proteinase, and levels of amylase and cellulase were in the order of ZRZ_ V > ZRZ_ VI > ZRZ_III (Supplementary Tables 1-3); the genes exhibited the same trend (Fig. 10). There were also different levels of esterase and chitinase, which are thought to be closely related to changes in the feeding habits of the larval B. areolata. Agrawal et al. 75 compared the activities of proteases, amylases, and lipases in carnivorous fish Wallago, omnivorous fish Clarias, and herbivorous fish Labeo, and found that the activity of carbohydrases is higher in the Labeo (herbivorous fish) than in the Wallago (carnivorous fish) and Clarias (omnivorous fish). Regarding proteases, maximum activity was found in Wallago (carnivorous fish). The difference in lipase activity was not as pronounced. Proteases and lipases exhibited the same activity in the order of carnivorous > omnivorous > herbivorous. Conversely, amylase activity was in the order of herbivorous > omnivorous > carnivorous. These results showed that there is a correlation between the normal diet of fish and the relative activity of digestive enzymes 75 . These findings are similar to the results of this study, which indicated that the activities of different digestive enzymes correlated with different feeding habits.
Immune-related protein at the metamorphosis stage. Mollusks lack an adaptive immune system, and mainly depend on the innate immune response as a defense against various pathogens. In this study, a variety of immune proteins were identified including C-type lectin, lysozyme, stress protein A, GST, Hsp70, cathepsin L and marginal zone and B1 B cell-specific protein-like (as part of the natural immune memory). C-type lectin, GST, and cathepsin L levels were in the order of ZRZ_III > ZRZ_V > ZRZ_VI, and lysozyme was highly expressed at ZRZ_V. The abundance of stress protein A was in the order of ZRZ_VI > ZRZ_V > ZRZ_III. Hsp70 was only detected at ZRZ_VI (Supplementary Tables 1-3). Regarding the mRNA expression level, cathepsin L was upregulated from ZRZ_III to ZRZ_VI, the inconsistency at the mRNA and protein levels provokes that investigations are needed to reveal the changes of cathepsin L genes during metamorphosis. GST was downregulated from ZRZ_III to ZRZ_VI. The expression levels of C-type lectin, lysozyme, stress protein A, Hsp70, marginal zone and B1 B cell-specific protein-like genes were not statistically different among the three stages (p > 0.05), C-type lectin and stress protein A were slightly upregulated at the ZRZ_V stage. Metamorphosis is the critical stage when larvae increase the expression of immune-related genes and respond to environmental stimulus 1 . Balseiro et al. 1 showed that at this metamorphic stage in Mytilus galloprovincialis, the expression levels of immune-related genes were higher than those in oocytes, suggesting that immune-related genes are active in mussel larvae. The upregulation of innate immune-related genes during metamorphosis was previously reported in ascidian Boltenia villosa 76 . Protein expression began to increase after the veliger stage to prepare larvae for settlement. The activation of innate immunity genes during metamorphosis may not only reflect maturation of the innate immune system but also might correlate with the resorption and re-structuring of larval tissues, as well as the ability of larval to detect and respond to bacterial settlement cues 1,76 . C-type lectins play important roles in the immunity of invertebrate. Intracellular lectins mainly play roles in protein trafficking and sorting, whereas extracellular lectins are The authors found that the expression level increased in the gastrulae, trochophore, and early D-shaped larvae (veliger stage), and then decreased in D-shaped larvae, and increased hundreds of times in metamorphosing larvae. In another study, lysozymes were selected for expression profile analysis throughout the different developmental stages (trochophore, veliger, metamorphosis, post-settlement, and spat) in mussel M. galloprovincialis, and it was observed that expression of lysozyme genes increased during the transition from trochophore to spat 1 . In this study, C-type lectin and lysozyme were highly expressed at ZRZ_III and ZRZ_VI, respectively, similar to previous studies. In addition, Hsps have dual functions, and are involved in the stress response and developmental processes 79 . Gunter et al. 79 observed higher levels of Hsp70 expression in tissues undergoing larval morphogenesis in the marine gastropod H. asinina. Hsp70 is expressed in unique and overlapping patterns in the prototroch, foot, and mantle, and returns to lower levels after morphogenesis is complete. These patterns of Hsp70 expression in H. asinina are similar to those observed in ecdysozoans and deuterostomes 79 . In this study, Hsp70 was only detected at ZRZ VI, akin to previous research studies.
Proteins related to attachment and metamorphosis. In marine invertebrates, there are usually several metamorphosis events in early embryo development, which include development of the larvae internal structure, external morphology, physiological processes, and behavior. In planktonic life at the early veliger stage to the benthic life with gradually disappearing velum, the larvae undergo the attachment process. The nervous system plays an important role in the activation and regulation of the larval metamorphosis. GABA is a neurotransmitter that has been widely used in the chemical induction of marine shellfish larvae metamorphosis 80 . Fang et al. 81 proposed in their studies on larval Arca granosa L. that tryptophan, GABA, and acetyl choline within certain concentration ranges exert clear inducing effects. GABA may significantly induce the larval attachment of H. discus 82 , which is consistent with the findings of Morse et al. 83 on H. rufescens. GABA induced both settlement and metamorphosis in the mussel M. galloprovincialis, the clams Venerupis pullastra and Ruditapes philippinarum, and the oyster Ostrea edulis 84 . In this study, the GABA-ergic synaptic signaling pathway was also significantly enriched (Fig. 6), which was likely related to larval attachment. The immune-related genes may be involved in the interaction with biofilm during settlement (i.e., attachment to a surface and metamorphosis into juveniles) when the larvae are searching for the best habitat to settle. Some studies have discussed the processes of metamorphosis in some marine organisms, and lectins (sugar binding proteins and glycoproteins) have been suggested to be involved in settlement, metamorphosis, and tissue remodeling 76,85-87 . Maki and Mitchell 88 showed that lectins are involved in the settlement and metamorphosis of marine invertebrate larvae, and proposed a biochemical lectin model for settlement and metamorphosis of some marine invertebrate larvae. The authors demonstrated that lectin on the surface of larvae "recognizes" and binds to a glycoconjugate in the exopolymer of marine bacteria. Matsutani et al. 89 showed that lectin-like factors are involved in larval settlement and metamorphosis in the abalone H. discus hannai. Bao et al. 78 Table 6. Comparison of mRNA and protein expression levels for 17 genes. Note: The letter "a" showed P > 0.05 between three development stages; the letter "b" showed P > 0.05 between ZRZ_III and ZRZ_V; the symbol "*" showed different variations between proteins and mRNA expression pattern.
SCiEntifiC REPORTS | (2018) 8:6269 | DOI:10.1038/s41598-018-24645-z been reported that calmodulin participates in shell formation 90 . Furthermore, a large amount of proteins related to calcium metabolism may be involved in calcium deposition. Regarding mRNA expression level, EF-hand calcium-binding domain-containing protein gene were upregulated at the ZRZ_III stage, and calcium-binding mitochondrial carrier protein SCaMC and calmodulin were slightly upregulated at the ZRZ_III and ZRZ-V stage respectively, suggesting that these genes were involved in shell formation. Whereas, in real-time PCR assay, the gene of EF-hand calcium-binding domain-containing protein had different patterns. Thus, it still remains unclear whether transcription/translation inhibitors would affect the metamorphosis of B. areolata.
Proteins related to body torsion and neural development. In this study, 14-3-3 and SCO-spondin were identified. Comparison of the ZRZ-III and ZRZ-V stages showed that 14-3-3 was significantly upregulated at ZRZ-V ( Supplementary Tables 1-3). 14-3-3 participates in growth and development, tumorigenesis and receptor signal regulation, molecular chaperone, and activation or inhibition of the target protein activity [91][92][93] . The overexpression of 14-3-3 in the body plays an important role in nerve cell differentiation 94 . Studies on brain development in humans and mice showed that 14-3-3 deficiency may cause brain death or severe brain deformity 95 . In the embryonic development of frogs, 14-3-3 in the early development stage participates in determining the left and right in the axial direction of the embryo 96 . The body torsion of the larval B. areolata is likely related with 14-3-3, and 14-3-3ε is related with the neural development of B. areolata. Additionally, SCO-spondin is required for neurogenesis during early brain development, Vera et al. 97 suggested that SCO-spondin is a vital embryonic cerebrospinal diffusible factor regulating the balance between proliferation and differentiation of the brain neuroepithelial cells. For 14-3-3 protein and SCO-spondin, the inconsistency at the mRNA and protein levels provokes the question of whether two gene participates in body torsion and neural development of B. areolata, that still suggest a further investigation.
Other proteins. In this study, the radial spoke head protein was identified, which is structural components of the axoneme and are suggested to be involved in producing flagellar patterns through regulation of the dynein arms 28 . The radial spoke head protein decreased from ZRZ_III to ZRZ_VI, indicating that it might be related to velum atrophy. Kyphoscoliosis peptidase plays a vital role in muscle growth 98 , and at the protein and mRNA level, the pattern was ZRZ_V > ZRZ_III > ZRZ_VI, indicating that it might correlate with muscle growth. We also noticed a discrepancy between the protein level and mRNA expression of the EF-hand calcium-binding domain-containing protein, cathepsin L, 14-3-3 protein, and SCO-spondin (

Conclusions
This study provides the first comparative analysis of proteomic profiles at different embryonic development stages of B. areolata. A total of 5583 proteins were identified. The DEPs were annotated and enriched in the GO and KEGG databases. Identification of these proteins will provide a basis for exploring larval development and metamorphosis mechanisms. We identified the protein profiles of B. areolata development related to shell formation, body torsion, changes in feeding habits, attachment and metamorphosis, and immune response.

Materials and Methods
Sampling. B. areolata were purchased from the Dongshan farm, Fujian province, selected individuals with vitality. Snails were maintained in an ecological flow-through water pond (20.0-30.0 °C, salinity 28.3-35.1, pH 7.8-8.5, oxygenated). The snails were fed daily with oysters and chopped fresh fish. Embryos were cultured, and fertilization, hatching, and developmental conditions were measured. Three groups of larval B. areolata at different developmental stages including the middle veliger stage before attachment (ZRZ-III), later veliger stage (velum atrophy) (ZRZ-V), and juvenile stage (ZRZ-VI) were used. Each group was replicated three times. Early veliger larvas were fed daily with algae. The late veliger stage sinks to the substratum and leads a benthic life, the feeding habits also gradually changed from phytophagy to sarcophagy, the snails were fed daily with oysters. Larval B. areolata were sampled from a culture pond containing 2,000,000-3,000,000 larvae. The sampled larvae were flushed with clean seawater, transferred to a 1.5 mL centrifuge tube with a suction pipe, rinsed twice with normal saline, frozen in liquid nitrogen, and stored at −80 °C.

Extraction and Quantification of Proteins.
Samples were thawed at 4 °C. Subsequently, 60 mg sample was weighed. Then 200 μL SDT lysis buffer (4% SDS, 100 mM Tris-HCl, 1 mM DTT, pH7.6) was added to the sample, followed by solubilization for 20 s with a tissue homogenizer; this was repeated three times. Then, 500 μL SDT lysis buffer is added. The mixture was shaken and heated for 10 min in a water bath. After 30 min ultrasonic treatment, the mixture was centrifuged for 30 min at 12,000 × g and 4 °C, and the supernatant was saved for subsequent experiments. The BCA assay was used for protein quantification.

SDS-PAGE for Protein Separation.
Approximately 20 μg of each protein sample was resolved and mixed with 5× loading buffer added at 1:5 ratio. The mixture was boiled for 5 min in a water bath and centrifuged for 10 min at 14,000 × g. The supernatant was used, and the mixture was electrophoresed for 60 min on 12.5% SDS-PAGE gels at 15 mA. Protein bands were visualized by Coomassie Blue R-250 staining. The mixture was transferred to a 30 kDa centrifuge tube for ultrafiltration and centrifuged for 15 min at 14,000 × g, after which the filtrate was discarded. After adding 200 μL UA buffer, the mixture was centrifuged for 15 min at 14,000 × g, and the filtrate was discarded. Then IAA (50 mM IAA in UA) was added in UA buffer (100 μL), and the mixture was oscillated for 1 min at 600 rpm. The mixture was allowed to settle for 30 min in the dark at room temperature and centrifuged for 10 min at 14,000 × g. UA buffer (100 μL) was added, and the mixture was centrifuged for 10 min at 14,000 × g; the procedures were repeated twice. Next, 100 μL of NH 4 HCO 3 buffer was added, and the mixture was centrifuged for 10 min at 14,000 × g, the procedures were repeated twice. Subsequently, 40 μL trypsin buffer (2 μg Trypsin in 40 μL NH 4 HCO 3 buffer) was added. The mixture was oscillated for 1 min at 600 rpm and allowed to settle for 16-18 h at 37 °C. The mixture was transferred into a new collecting tube and centrifuged for 10 min at 14,000 × g, after which the filtrate was collected. The filtrate was desalinized with C18 SD Extraction Disk Cartridge and quantified at OD280.

LC-MS/MS Analysis on the Protein Digestion Product.
Approximately 2 μL product of protein digestion was used for LC-MS/MS analysis and separated using a nanoliter HPLC EASY-nLC1000 system. The mobile phase A was 0.1% formic acid solution with 2% acetonitrile, and the mobile phase B was 0.1% formic acid solution with 84% acetonitrile. The chromatographic column EASY column SC200 150 μm * 100 mm (RP-C18) was balanced with 100% A solution. The sample was loaded onto the EASY column SC001 traps 150 μm * 20 mm (RP-C18) through an automatic sampler and separated by a chromatographic column at the flow rate of 400 nL/ min. The gradient elution procedure was as follows: the percentage change of the mobile phase B for 0-100, 100-108, and 108-120 min is 0-45%, 45-100%, and 100% of the linear change, respectively. After the capillary separation by HPLC, the product of protein digestion is determined using Q-Exactive mass spectrometer with the following parameters: duration: 120 min; detection method: positive ion detection; parent ion scanning scope: The samples from each of the three stages, namely ZRZ-III, ZRZ-V, and ZRZ-VI, were prepared in triplicate. Analysis was performed by one-way ANOVA, and P < 0.05 were considered statistically significant.

GO Functional Annotations.
The steps for GO annotation were as follows 99 . Target proteins were assembled and prepared through sequence alignment and batch extraction of target protein and BLAST of GO terms (mapping). GO annotation is a process of scoring the similarity between the target sequence and the sequence to be aligned to determine the reliability of the source of the GO term and the GO term structure in directed acyclic graphs. Only the GO term that satisfies the preset score may be annotated to the target protein sequence. For the supplementary annotation augmentation, the conservative motif in the EBI database of InterProScan is aligned with the target protein. The motif-related functional information is annotated to the target protein sequence. The annotation is further supplemented with ANNEX. The links between different GO terms were established, thereby improving the accuracy of the annotation. Consequently, the GO annotation results are obtained.

GO Enrichment Analysis on Deps.
Fisher's exact test was used to evaluate the significance level of protein enrichment under each GO term.
KEGG Pathway Annotations. The steps were as follows. The protein sequence is aligned with the KEGG database 100 , and homologous KEGG genes are obtained. The resulting gene was screened through Bi-Directional hit rate, and the orthologous candidate genes were obtained. The KO grouped data were obtained by KO based on the probability and Heuristic scoring, and the KO ranking chart is obtained. Finally, KEGG pathway mapping were finished. KEGG Pathway Enrichment Analysis on Deps. The KEGG pathway enrichment analysis was similar to the GO enrichment analysis. In this analysis, the KEGG pathway was the unit, and all of the qualitative proteins were background. The significance level of protein enrichment of each pathway was calculated by Fisher's exact test to determine the significantly affected metabolism and the signal transduction mode.
Experimental validation using qPCR. The differentially expressed protein genes were validated by qPCR to confirm the proteomic results, which was conducted using a QuantStudio 6 Flex Real-Time PCR System with Thermo Scientific DyNAmo ColorFlash SYBR Green qPCR Kit according to manufacturer instructions. Primer sequences were designed based on each identified gene sequence using Primer Premier 6 software (Premier Biosoft, USA) (Supplementary Table 4). PCR amplification experiments were conducted in triplicate SCiEntifiC REPORTS | (2018) 8:6269 | DOI:10.1038/s41598-018-24645-z under the following conditions: Initial denaturation, 95 °C for 7 min, then denaturation 40 cycles of 95 °C for 10 s, Annealing/extension 60 °C for 30 s. The results were normalized using actin (c264956_g1; gi|404452331|g-b|AFR74960.1| actin, partial [Rana clamitans]) and EF1-F (c233532_g1; gi|524916452|ref|XP_005113003.1| PREDICTED:elongation factor 1-alpha [Aplysia californica]) for each sample and the 2 −ΔΔCT method. The expression levels of each gene in the three developmental stages were compared using a two-sided Student's t-test, and differences were considered statistically significant at p < 0.05.