Proteome of larval metamorphosis induced by epinephrine in the Fujian oyster Crassostrea angulata

Background The Fujian oyster Crassostrea angulata is an economically important species that has typical settlement and metamorphosis stages. The development of the oyster involves complex morphological and physiological changes, the molecular mechanisms of which are as yet unclear. Results In this study, changes in proteins were investigated during larval settlement and metamorphosis of Crassostrea angulata using epinephrine induction. Protein abundance and identity were characterized using label-free quantitative proteomics, tandem mass spectrometry (MS/ MS), and Mascot methods. The results showed that more than 50% (764 out of 1471) of the quantified proteins were characterized as differentially expressed. Notably, more than two-thirds of the differentially expressed proteins were down-regulated in epinephrine-induced larvae. The results showed that “metabolic process” was closely related to the development of settlement and metamorphosis; 5 × 10− 4 M epinephrine induced direct metamorphosis of larvae and was non-toxic. Calmodulin and MAPK pathways were involved in the regulation of settlement of the oyster. Expression levels of immune-related proteins increased during metamorphosis. Hepatic lectin-like proteins, cadherins, calmodulin, calreticulin, and cytoskeletal proteins were involved in metamorphosis. The nervous system may be remodeled in larval metamorphosis induced by epinephrine. Expression levels of proteins that were enriched in the epinephrine signaling pathway may reflect the developmental stage of the larvae, that may reflect whether or not larvae were directly involved in metamorphosis when the larvae were treated with epinephrine. Conclusion The study provides insight into proteins that function in energy metabolism, immune responses, settlement and metamorphosis, and shell formation in C. angulata. The results contribute valuable information for further research on larval settlement and metamorphosis. Graphical abstract


Background
Oysters are a group of commercially important species cultured along the coast of China. The annual production of oysters was about 4.88 million tonnes in 2017 (China Fishery Statistical Yearbook for 2018) [1]. The Fujian oyster Crassostrea angulata accounts for about 50% of the total oyster production in China [2,3].
Embryonic and larval development are key phylogenetic events [4], and embryological research contributes to the development of environmental pollution monitoring and the aquaculture industry [5]. Size and age at reproduction of the offspring may have important consequences for population dynamics and demography [6]; larval developmental plasticity is a crucial source of variation and can directly influence the evolution of populations and species in adult phenotypes [7,8]. Oysters are a typical two-phase life cycle species that have a pelagic phase and a benthonic phase; the transition between the two forms usually occurs rapidly [9]. Larval settlement and metamorphosis are importmant transition periods associated with the evolution of mollusks, and these factors also impact population distribution, phenotypic differentiation, and speciation [10,11]. Even in cultured oysters, metamorphosis is a crucial step for the aquaculture facility, because oysters that have successfully settled usually show increased survival.
Metamorphosis of oysters includes complex morphological and behavioral changes that are irreversible and that are essential for survival. Pre-settlement or postmetamorphosis larvae are proposed to be regulated by endogenous chemical signals and endocrine proteins. However, the molecular mechanisms underlying metamorphosis of model species, such as frogs and fruit flies, have been studied in much greater depth than those in marine mollusks. In the Fujian oyster, Cacaspase-2, Cacaspase-3 mRNA, and settlement and metamorphosis-related protein (SMRP1) are highly expressed in larvae during settlement and metamorphosis [12,13]. Cranfield et al. described the behaviour of Ostrea edulis before attachment [14][15][16]. Ke and Feng [17] described the settlement process of Crassostrea angulata. The choice of settlement substratum is critical because an inappropriate settlement substratum is often fatal to juveniles [18]. Scallop or oyster shells are often used as substrata in oyster farms; the shells are placed into cages for raft-or long-line culture [19]. In addition, the traditional farming methods are laborous and costly; food acquisition and growth rate of the oysters are often affected by the high-density farming that causes large individual differences and reduced economic value. Thus, it is important to understand the breeding biology of oysters.
Oyster hatcheries have been developing multiple methods for clutchless oyster spat production for many years [20]. The most widely used inducer for oysters is epinephrine (EPI, also known as adrenaline) [21]. Epinephrine induction has shown the capability to generate non-attached spat (juveniles) in various oyster species [20,22,23]. Coon et al. [20] report that EPI and norepinephrine are able to induce cultchless spats, demonstrating that EPI is able to induce oyster metamorphosis without settlement or negative effects on survival and development [20]. In comparison with oysters attached to the substratum, clutchless oysters are superior in shape and are more uniform, traits that are economically favorable [24]. Thus, the use of EPI has become a common practice in many bivalve hatcheries [25,26]. EPI has been widely used to induce the settlement and metamorphosis of bivalve mollusks [25,[27][28][29][30]. The positive effects of EPI on metamorphosis in a wide range of bivalve species have been summarised recently [31].
Nonetheless, the molecular mechanisms of EPI induction are unclear. Coon et al. proposed that adrenergic receptors are involved in the regulation of settlement and metamorphosis in oysters [20,32]. In 2012, molecular characterization and functional analysis of adrenergic-like receptors during larval metamorphosis in Crassostrea angulata was studied by our research group [33].
Recently, many studies of the molecular mechanisms of metamorphosis have used sequencing techniques, transcriptome, and genomic data analysis [34]. However, genetic data are unable to reveal the complete molecular mechanisms [35], because there is no direct correlation between gene expression intensity and protein abundance [36]. Proteomic analysis is a crucial approach for studies of bivalve oyster metamorphosis [35].
Proteomics analysis has been carried out on multiple marine genera during their settlement and metamorphosis stages [37][38][39][40][41][42]. As such, proteomic analysis has become a powerful tool to investigate the molecular mechanisms involved during different stages of development of marine invertebrates, and the method has been employed in embryonic developmental stages of several mollusks, including larvae [35,36,43,44] reproduction [45], attachment, and metamorphosis [46][47][48][49][50]. These studies provide comprehensive information on the molecular mechanisms involved. Proteomic analysis is an appropriate approach to study larval settlement and metamorphosis of C. angulata.
Label-free quantification proteomics involves the application of a label-free quantification algorithm used in Maxquant for proteomics data. This technique is rapid, clean, and simple. The method can directly and precisely quantify protein expression without using labeling, and it has been broadly applied to biomarker discovery and proteomic profiling [51]. In this study, label-free proteomic analysis was carried out on larvae of the oyster C. angulata that had been induced by epinephrine during settlement and metamorphosis. The collected data were analysed using Scaffold Proteomics Software with database [52]. In addition, gene function annotation, GO enrichment, KEGG enrichment, and transcriptome and proteome combined analysis and qPCR verification were carried out. The results provide insights into the different protein groupings active during energy metabolism, shell formation, attachment and metamorphosis, and immune-related activities in C. angulata.

Identification and quantitative results
In this study, a total of 412,946 spectra were identified (Supplementary Fig. 1), with 137,372 unique specta, 102, 696 peptide fragments, 60,245 unique peptide fragments, and 1471 proteins. The relative molecular weights of the identified proteins were mainly distributed in the range of 10-60 kDa, of which 20-30 kDa was the largest group (Supplementary Fig. 2A), accounting for 17.3%; Supplementary Fig. 2B showed that the length of the identified peptide segments was mainly distributed in 10-18 kDa, of which 11-14 kDa was the largest group, accounting for 35.0%. Supplementary Fig. 2C illustrates the coverage of the identified peptide segments, in which 49.4% coverage was more than 20%. Supplementary Fig. 2D shows that a number of identified proteins contained only one or two peptide segments, and three peptide segments accounted for 48.7%.

Identification and analysis of differentially expressed proteins (DEPs)
Protein abundances with at least 1.5-fold change among the groups were compared and considered as upregulated or down-regulated (Fig. 1a). A total of 611 proteins were differentially expressed in at least one larval stage among the PL-PA-MET process. There were 117 co-expressed proteins in the PL-PA and PA-MET groups ( Fig. 1 b). The results of a cluster analysis showed that expression levels of most proteins had a consistent trend in the PL-PA-MET developmental process ( Fig. 1 b).
There were 246, 341, and 230 DEPs in the PL-eSEN group, the PL-einSEN group, and the eSEN-einSEN groups, respectively (Fig. 1a). There were 35 DEPs present in all 3 groups, of which 22 were down-regulated in the PL-eSEN and PL-einSEN groups, and 3 proteins were down-regulated in the eSEN group (Fig. 2a) and upregulated in the einSEN group (Cluster 1 of Fig. 2c). A total of 101 proteins were differentially expressed only in the PL-einSEN group and eSEN-einSEN group ( Fig. 2a-b); there were 26 proteins that were down-regulated in the eSEN group but up-regulated in the einSEN group. The trend of expression is shown in cluster 2 of Fig. 2c.
There were 560, 539, and 443 DEPs in the PL-MET, PL-eMET and MET-eMET groups, respectively, of which 286, 235, and 194 proteins were up-regulated and 274, 304, and 249 proteins were down-regulated, respectively (Fig. 1a). There were 97 DEPs common to all three groups (Fig. 3a). Among these, nine proteins were down-regulated in MET and up-regulated in eMET. The trend of expression is shown in cluster 3 of Fig. 3c, while 13 proteins were up-regulated in MET and downregulated in eMET, as shown in cluster 4 of Fig. 3c. There were 132 DEPs observed only in the PL-MET and MET-eMET groups, and 116 DEPs found only in the PL-eMET and MET-eMET group (Fig. 3a). In addition, 197 proteins were differentially expressed in the PL-MET and PL-eMET groups but not in the MET-eMET group. Among these, 118 of the up-regulated proteins may be related to metamorphosis. Figure 4a displays the Venn diagram of DEPs in five groups, with the PL group as a control group. A total of 954 proteins were differentially expressed in at least one group. Interestingly, 54 proteins were differentially expressed in all five of the groups. The hierarchical cluster analysis (Fig. 4b) showed that the expression trends of most proteins were basically similar. In addition, the PL-MET and PL-eMET groups had the highest numbers of specific differentially expressed proteins, 159 and 126 proteins, respectively, and the PL-eSEN group had the least number, only 21 differentially expressed proteins. Table 1 shows some DEPs, including long chain fatty acid CoA ligase (ACSBG), calnexin, and glutathione Stransferase.     , mitogenactivated protein kinase 1-like, soma ferritin-like, universal stress protein A-like protein were higher expression in eMET group. Meanwhile, calcium uniporter protein, calcium-binding mitochondrial carrier protein Aralar1, calcium-transporting ATPase, calmodulin, calnexin, calumenin-like isoform X1, cathepsin L1-like, collagen alpha, drebrin-like protein B isoform X2, laminin subunit alpha (GI405969732,762,109,423), laminin subunit gamma-1, vinculin-like isoform X7, lethal (2) giant larvae protein (GI 762141840, 405,952,027), neurexin-4, neuroglian, neuroglian-like isoform X1, neutral ceramidase-like, HEAT repeat-containing protein 7A, heat shock 70 kDa protein 12B, heat shock 70 kDa protein 14-like, heat shock protein 27-like, Ras-like GTPbinding protein RHO, ras-like protein 3 isoform X2, Rab-7a, Rho GTPase-activating protein 17, V-type proton ATPase subunit H-like isoform X3 were lower expression in eMETgroup (Supplementary Table 4).
The differentially expressed proteins with important physiological functions among six groups along with calculated protein volumes are shown in Supplementary  Table 5. For statistical analysis, the expression levels of proteins with the same name were merged. The results are shown in Supplementary Table 6.

Energy metabolism proteins
In present study, we identified several energy metabolism proteins. Among them, ATP synthase, D-3phosphoglycerate dehydrogenase-like, enolase and malate dehydrogenase were not statistically different among the six groups. 6-phosphogluconate dehydrogenase, decarboxylating, glucose-6-phosphate 1-dehydrogenase, glucose-6-phosphate isomerase-like, succinate dehydrogenase, pyruvate dehydrogenase E1 component subunit alpha type II, mannose-6-phosphate isomerase-like, citrate synthase and mitochondrial-like isoform X1 were up-regulated in PL group and down-regulated in MET. Pyridoxal-dependent decarboxylase domain-containing protein 1 was up-regulated in MET. Most proteins involved in carbohydrate metabolic processes were upregulated in PL group and down-regulated in MET group. (Supplementary Table 6).

Immune proteins
In this study, a variety of immune proteins were identified. Among them, cathepsin Z, cathepsin L, Hsp60, Hsp70, Hsp27 and IgGFc-binding protein-like were upregulated in MET, these proteins were also highly expressed in eMET. Hepatic lectin-like and Cu/Zn-SOD -like isoform X1 were up-regulated in eMET, these proteins were also highly expressed in MET. On the contrary, galectin-6, glutathione-S-transferase, universal stress protein A-like protein, galectin-9-like isoform X2 and Mn-SOD were down-regulated in eMET; universal stress protein A-like protein, galectin-9-like isoform X2 and Mn-SOD also were down-regulated in MET (Supplementary Table 6).

Shell formation proteins
In this study, cadherin, calmodulin, calreticulin, caltractin-like and calumenin-like isoform X2 were upregulated in MET and eMET groups; CaM kinase was down-regulated in MET and eMET groups; calcium uniporter protein was down-regulated in eMET group; Calnexin was up-regulated in MET, down-regulated in eMET group; EF-hand calcium-binding domaincontaining protein was down-regulated in MET group, up-regulated in eMET group; Troponin T was upregulated in eMET group. The diversity of expression patterns also indicates the complexity of shell formation regulation.

GO annotation results
All of the identified proteins were annotated by GO. A total of 11,028 proteins were annotated, accounting for 55.1% of the total identified proteins. Among these, 55 secondary GO terms were classified ( Supplementary Fig. 3). From a total of 51,618 nodes, 18,814 were involved in biological processes, mainly in cell processes and metabolic processes, with 5405 and 5247 nodes accounting for 28.7 and 27.9% of the total GO term nodes, respectively. For cellular components, there were 20,849 nodes, mainly distributed in cell and cell part classifications. In the molecular function classification, there were 11,955 nodes, mainly involving two types of GO terms: binding and catalytic activity, with 5718 and 4547 nodes, accounting for 47.8 and 38.0%, respectively.

KEGG annotation results
The identified proteins were annotated to KO according to NR Gene ID. A total of 5742 proteins were annotated, accounting for 28.7%. In primary pathway classification, the number of proteins participating in metabolism, genetic information processing, environmental information processing, cellular processes and organismal system were 1258, 630, 996, 660, and 1529 proteins, respectively (Supplementary Fig. 4), accounting for 24.8, 12.4, 19.6, 13, and 30.1%, respectively. There were 32 secondary KO classifications ( Supplementary Fig. 4). The pathway annotated to the signal transduction was the largest group, with 877 proteins accounting for 17.3%, followed by endocrine system with 435 protein annotations, accounting for 8.6%. The numbers of proteins participating in immune system, carbohydrate metabolism, nervous system and amino acid metabolism were 289, 279, 278, and 269 proteins, accounting for 57.0, 55.0, 54.8, and 53.0% of the total annotated proteins, respectively.

GO enrichment analysis of DEPs
For GO significance analysis of DEPs, a P-value < 0.05 and FDR < 0.05 were used as a threshold to select significant GO terms. Table 2 shows the GO enrichment in all of the groups; Supplementary Table 7 shows the extremely significant GO enrichment in the PL-PA and PA-MET groups. In PL-PA and PA-MET groups, 21 GO terms were extremely significantly enriched, including 14 GO terms in biological processes, 4 GO terms in cell components, and 3 GO terms in molecular functions (Supplementary Table 7).
The statistical analysis of GO enrichment in the PL-eSEN, PL-einSEN, and eSEN-einSEN groups is showed in Table 2. The number of DEPs with significantly enriched GO terms in the PL-einSEN group was the largest, while the eSEN-einSEN group had the least number of DEPs and the least significantly enriched GO terms. Supplementary Table 8 shows that seven extremely significant GO terms were only expressed in the PL-eSEN and eSEN-einSEN groups; these were regulation of cell adhesion, histone lysine methylation, F-actin capping protein complex, neuropeptide hormone activity, dimethylallyltranstransferase activity, chaperone binding, and peptide-methionine (R)-S-oxide reductase activity. GO entries that were extremely expressed only in the PL-einSEN and eSEN-einSEN groups are shown in Supplementary Table 8.
The statistical analysis of GO enrichment in PL-MET, PL-eMET, and MET-eMET groups is shown in Table 2.
The significantly enriched and extremely significantly enriched GO terms were involved in the molecular process among three groups; the most extremely significant GO terms were in the PL-MET group, which had 37 terms. Supplementary Table 9 shows that GO terms that were extremely significant only in the PL-MET and MET-eMET groups, only in the PL-eMET and MET-eMET groups, and only in the PL-MET and PL-eMET groups.
In addition, with PL as a control group, precorrin-2 dehydrogenase activity was the only extremely significantly enriched GO term for molecular functions in all five of the groups.
KEGG enrichment analysis of DEPs P < 0.05 was used as a determiner of significant enrichment in the KEGG pathway analysis. Pathways were classified in the PL-PA group, the PA-MET group, the PL-eSEN group, the PL-einSEN group, the eSEN-eMET   Table 10). The pathways related to carbohydrate metabolism included propanoate metabolism, citrate cycle (TCA cycle), butanoate metabolism, and pentose and glucuronate interconversion. The pathways related to lipid metabolism included fatty acid elongation, sphingolipid metabolism, biosynthesis of unsaturated fatty acids, fatty acid degradation, and primary bile acid biosynthesis (Supplementary Table 10). In addition, there were eight significantly enrichment pathways in both PL-PA and PL-MET groups.
From the pediveliger of the oyster to post-epinephrineinduced metamorphosis, in the PL-eSEN group, 14 DEPs were enriched in 17 pathways, of which 7 were related to organismal systems, and 4 were related to metabolism (Fig. 5). Compared with PL-eSEN group, 15 of 32 DEPs in the eSEN-eMET group were related to organismal systems, and 12 were related to metabolism. There were 12 unique pathways related to organismal systems in the eSEN-eMET group (Supplementary Table 11), including adrenergic signaling in cardiomyocytes, adipocytokine signaling pathway, and thyroid hormone synthesis; seven pathways were related to metabolism (Supplementary  Table 11), including glycosphingolipid biosynthesisganglio series, alanine, and aspartate and glutamate metabolism. In addition, there were six pathways with significant enrichment in both the PL-eSEN and eSEN-eMET groups.
Twenty-three DEPs were significantly enriched in 50 pathways in the PL-einSEN group. Twenty-eight were related to organismal systems pathways (Fig. 5). Twenty-five were unique to the PL-einSEN group (Supplementary Table 12). Many of the DEPs were related to hormones or neurotransmitters, including insulin signaling pathway, adrenergic signaling in cardiomyocytes, GnRH signaling pathway, cholinergic synapse, dopaminergic synapse, and other pathways (Supplementary Table 12). The identified proteins involved in adrenergic signaling in cardiomyocytes are shown in Table 3. In the MET-eMET group, 26 DEPs were enriched in 38 pathways, 7 of which were related to organismal systems pathways. Twenty-four were related to metabolic pathways, of which 10 were related to carbohydrates and lipids, accounting for 41.7% (Fig. 5).

mRNA expression of six differentially expressed genes
To investigate the transcriptional levels of identified proteins, six DEPs were chosen for qPCR verification (Fig. 6) (Calm-calmodulin-like; PSMB-proteasome subunit beta; HYOU-hypoxia up-regulated 1; GST-glutathione Stransferase omega; PLCPI-Leukocyte cysteine proteinase inhibitor 1; and PTPN11-tyrosine-protein phosphatase non-receptor type 11). Compared with the protein expression patterns, the mRNA expression of two of the six genes (Calm and PTPN11) exhibited similar patterns; three genes (PSMB, HYOU, and GST) were differentially expressed at the transcriptional level but showed no significant changes in the protein levels. Figure 7 was obtained from the http://string.embl.de/ website and shows the predicted interactions of the identified DEPs. Major clusters were associated with protein biosynthesis, motor protein, cilium movement, axoneme assembly, response to calcium ion, actinbinding, endocytosis, oxidoreductase, carbohydrate metabolism, mRNA processing, and DNA-binding (Supplementary Table 13).

Discussion
Studying the biological process of metamorphosis is the key to understanding the origin and evolution of mollusk development. The oyster transforms from a pelagic larva to a benthic adult, during which settlement and metamorphosis occur. Larvae undergo dramatic tissue transition [32], and the transition usually happens rapidly [9]. Some larval tissues degenerate during settlement and metamorphosis, e.g., the anterior adductor, the ciliated velum, and the ventral retractors [53]. The transition initiates histogenesis of some adult tissues (e.g., calcified shells and gill) as well as the formation of eyespot larvae (pediveliger)-specific tissues (e.g., eyespot and foot) during settlement and metamorphosis. As the oyster adult body plan is distinct from the larva, understanding how the molecular regulation is reorganized is crucial to gaining insight into the complexity of the transition [54]. In this study, a total of 1471 proteins were identified and quantified. Quantitative results showed that 764 DEPs were related to the larval developmental processes of settlement and metamorphosis, and 683 DEPs were related to larval development process of direct metamorphosis induced by epinephrine. Of the latter, there were 417 DEPs present in both the PL-PA-MET and PL-eSEN-eMET groups.

Comparison of eSEN (EPI-sensitive larvae) and einSEN (EPI -insensitive larvae) groups
There were 246 DEPs in the PL-eSEN group; the number of DEPs in the PL-einSEN group was nearly onethird more, at about 341 proteins. In this study, the larvae that were not successfully induced displayed an excited planktonic state after the epinephrine-inducing fluid was replaced, while the larvae that were successfully induced quietly sank to the bottom of the tank. This may be due to the different motion states between the PL-eSEN group and the PL-einSEN group, where the number of DEPs differed. In addition, downregulated proteins accounted for 67.1 and 66.6% of the differentially expressed proteins in the PL-eSEN group and PL-einSEN group, respectively. Most of the DEPs in this experiment were down-regulated. The reasons for this may have been as follows: firstly, the result may have been related to sampling. All of the larvae contracted their foot and velum rapidly when epinephrine was added, then closed their shells and sank to the bottom. At this stage, only the rotation of the visceral mass could be observed under the microscope, and most of the larvae remained in a dormant-like state for a long time. The results suggested that most of the physiological processes had slowed down or stopped. In this study, the epinephrine treatment time was only for 2 h. It is possible that larvae could not express new proteins necessary to complete the life cycle in such a short time.
The number of significantly enriched KEGG pathways in the PL-einSEN group was three times that in the PL-eSEN group. There were 50 significantly enriched pathways, including cytochrome P450-mediated substance metabolism, physiological rhythm transformation, extracellular matrix receptor function, and endoplasmic reticulum protein processing, present in both the PL-eSEN and PL-einSEN groups. Cytochrome P450 family was observed with a total of 491 genes belonging to P450 family members in pearl oyster [55]. Metabolism mediated by cytochrome P450 enzymes can decompose lipophilic toxic substances into aqueous solution [56], and this was up-regulated in eSEN and down-regulated in the einSEN group.
In addition, all of the proteins enriched in this functional module were down-regulated, including carbonyl reductase 1 (CBR 1), aldehyde dehydrogenase (E 1.2.1.5), and glutathione S-transferases (GSTs); the detoxification function of glutathione S-transferases has been widely studied [57]. In this study, compared with the PL group, the expression levels of GSTs were low in the einSEN and eSEN groups (Supplementary Table 6). This result suggested that epinephrine treatment for 2 h was not toxic to larvae, so proteins related to detoxification were down-regulated.
The unique signaling pathways in the PL-einSEN group included adrenergic signaling in cardiomyocytes, dopaminergic synapse, insulin signaling pathway, and  Table 12). The expression trend of most DEPs enriched in these pathways was down-regulation; there were eight DEPs in the adrenergic signaling in cardiomyocytes, among which six were down-regulated, including serine/threonine-protein phosphatase, cAMP-dependent protein kinase (PKA), mitogen-activated protein kinase (MAPK), calcium/calmodulin-dependent protein kinase II (CAMK 2), ryanodine receptor, and calmodulin (Table 3). Although the expression of these proteins was down-regulated in the eSEN group, only three proteins (calmodulin, ryanodine receptor, calcium/calmodulin-dependent protein kinase) Fig. 7 Predicted interactions of identified DEPs. Different line colours represent types of evidence for association. Major clusters associated with protein biosynthesis, motor protein, cilium movement, axoneme assembly, response to calcium ion, actin-binding, endocytosis, oxidoreductase, carbohydrate metabolism, mRNA processing, DNA-binding. Proteins without interactions have been removed from the graph. Protein abbreviations and corresponding full name are shown in supplementary Table 13 were differentially expressed (Table 3). We speculated that the expression levels of proteins related to the epinephrine pathway in the larvae of the eSEN group were higher than in the einSEN group before epinephrine treatment. After the same epinephrine stimulation, although the downregulation of differential proteins in the two groups was the same, the expression levels of proteins in the eSEN group were higher, so the einSEN group showed more DEPs. In this study, the proteins related to the epinephrine pathway were at higher expression levels in pediveligers (PL), and the expression of those proteins may reflect the developmental stage of the larvae, as this stage was directly involved in metamorphosis when the larvae were treated with epinephrine.

Proteins involved in energy metabolism in settlement and metamorphosis
There were 29 KEGG enrichment results related to metabolic pathways in the PL-MET group, including carbohydrate and lipid metabolism. The pathways related to carbohydrate metabolism included propanoate metabolism, citrate cycle (TCA cycle), butanoate metabolism, and pentose and glucuronate interconversions. The pathways related to lipid metabolism included fatty acid elongation, sphingolipid metabolism, biosynthesis of unsaturated fatty acids, fatty acid degradation, and primary bile acid biosynthesis (Supplementary Table 10). First, the metabolic processes can generate energy to ensure the normal development of oysters while also providing the material needed for oysters' development, as the processes of settlement and metamorphosis involve changes in organs and tissues. Carbohydrates and lipids provide energy and nutrients for embryonic development [58]. Isocitrate dehydrogenase is involved in energy homeostasis, and this was down-regulated in eMET (Supplementary Table 6). ATP synthase is needed for cellular energy interconversion. Enolase is a metalloenzyme responsible for the catalysis of the conversion of 2-phosphoglycerate to phosphoenolpyruvate in glycolysis. Glycolysis-related proteins play a key role in the embryonic development of zebrafish and ascidians [59]. Most proteins involved in carbohydrate metabolic processes were up-regulated in the PL group and down-regulated in the MET group (Supplementary Table 6). These results suggested that the metabolism of larvae was restored to a stable state after settlement and metamorphosis.
The KEGG enrichment of the MET-eMET group showed that KEGG pathways were mainly distributed in metabolic processes and organism systems. Lipid metabolism and polysaccharide synthesis and metabolism were the main metabolic pathways, including fatty acid elongation, fatty acid degradation, sphingomyelin metabolism, biosynthesis of unsaturated fatty acids, carbon metabolism, and glycosphingolipid biosynthesis. Jensen et al. [60] reported artificial induction of larval metamorphosis by free fatty acids in the polychaete Phragmatopoma californica, suggesting that free fatty acids induce larval metamorphosis by operating physiologically downstream or parallel to the natural inducer. The results indicated that the metamorphosis process of oysters requires significant energy consumption and that lipid metabolism is involved in larval metamorphosis.

Proteins involved in settlement between PL and PA groups
In marine invertebrates, the velum disappears gradually from planktonic life at the early veliger stage to the benthic form as the larvae undergo the attachment process. In 1996, Clare reported that a low concentration of Ca 2+ inhibited the settlement of Barnacle larvae [61]. Yamamoto et al. [62] and Zhang et al. [63] found that calmodulin (CaM) and calmodulin-determined myosin light chain kinase (MLCK) played important roles in the settlement process of barnacles. Zhang et al. [63] showed that CaM was highly expressed in cyprids in the swimming larvae that are competent to attach and undergo metamorphosis. In this study, in comparing the PL and PA groups, calcium/calmodulin-dependent protein kinase type II delta chain (GI 405964165) and calmodulin (GI20137620, GI762161385) were more highly expressed in the PL group. This result indicated that calmodulin was involved in the regulation of settlement of oysters.
Previous studies have found that the p38 MAPK pathway was activated during the settlement process of barnacles [64,65]. In this study, protein ID c100749_g1 was identified as MAP kinase-activated protein kinase 2-like (Supplementary Table 5). In comparisons of the PL and PA groups, mitogen-activated protein kinase 2-like was expressed at a higher level in the PL group. The NOTCH signaling pathway was involved in the regulation of somatic division of Capitella sp. [66]. In this study, the expression levels of neurogenic locus Notch protein in the PL, MET and eMET groups were 0.00, 4.77 and 13.89, respectively. The protein was up-regulated in the MET and eMET groups (Supplementary Table 6). This result indicated that MAPK and NOTCH signaling pathways were also involved in the regulation of settlement and metamorphosis of oysters. The classical signal pathway closely related to cell proliferation and apoptosis, the mitogen-activated protein kinase (MAPK) pathway, has also been well studied. Those pathways are involved in the regulation of settlement and metamorphosis of marine invertebrates [48,64,65,67].

The expression levels of immune-related proteins increased in metamorphosis
Metamorphosis is the key stage when the expression levels of immune-related genes increase and respond to environmental stimuli [68]. Balseiro et al. [68] revealed that the expression levels of immune-related genes were higher than those in oocytes at this metamorphic stage in Mytilus galloprovincialis. The results suggested that immune-related genes were active in mussel larvae. Previous study reported up-regulation of innate immunerelated genes during metamorphosis in the ascidian Boltenia villosa [69].
Due to the lack of adaptive immunity in mollusks, lectins act as determinants of phagocytosis [70]. Hsp70, as a member of the molecular chaperones, is a highly conserved stress protein that can protect cells from harmful assault. In addition, Hsps are involved in developmental processes [71]. Gunter et al. [71] found higher expression levels of Hsp70 in tissues undergoing larval morphogenesis in the marine gastropod H. asinine, similar to the patterns observed in ecdysozoans and deuterostomes [71]. Cu/Zn Superoxide dismutase (Cu/Zn-SOD) is closely related to immunity in mollusks. The enzyme can enhance the activity of phagocytic cells and the body's immune function and protect the cells from ROS poisoning [72]. In this study, cathepsin Z-like, cathepsin L, 60 kDa heat shock protein, heat shock 70 kDa protein 14-like (GI762129389), heat shock protein 27-like and IgGFc-binding protein-like were up-regulated in MET; subsequently, these proteins were highly expressed in eMET (Supplementary Table 6). Hepatic lectin-like (GI762109657) and superoxide dismutase (Cu-Zn)-like isoform X1 were up-regulated in eMET; subsequently they were highly expressed in MET (Supplementary Table 6). These results might indicate maturation of the innate immune system and the correlation between the resorption and restructuring of larval tissues and the ability of larvae to detect and respond to bacterial settlement cues [68,69]. Protein expression began to increase after the veliger stage in order to prepare the larvae for metamorphosis.

Proteins involved in metamorphosis
Some previous studies of metamorphosis in marine organisms have suggested that lectins are involved in settlement, metamorphosis and tissue remodeling [48,69,[73][74][75]. Maki and Mitchell [76] proposed a biochemical lectin model for settlement and metamorphosis of marine invertebrate larvae, and they demonstrated that lectin "recognizes" and binds to a glycoconjugate in the exopolymer of marine bacteria. Matsutani et al. [77] discussed how lectin-like factors were involved in larval settlement and metamorphosis in the abalone H. discus hannai. Bao et al. [78] suggested that a putative C-type lectin fold might play an important role during the metamorphosis of Mizuhopecten yessoensis. In this study, hepatic lectin-like (GI762109657) was up-regulated in the MET and eMET groups (Supplementary Table 6).
These results suggested that hepatic lectin-like proteins are involved in the metamorphosis of C. angulata.
In this study, cadherins, calmodulin and calreticulin were up-regulated in the MET and eMET groups (Supplementary Table 6). Cadherins are a family of homophilic Ca 2+ -dependent cell adhesion molecules controlling animal morphogenesis [79]. Chen et al. (2012) found that calmodulin affected settlement and metamorphosis of Hydroides elegans. In this study, the results were consistent with those from a study by Chen et al. on the expression of calmodulin in the processes of settlement and metamorphosis of the polychaete Hydroides elegans [80]. This suggests that calmodulin may play an important role in oyster metamorphosis.
Comparing the PA and MET groups, apoptosisinducing factor 3-like isoform X1 was expressed at a lower level in the MET group. Apoptosis has also been reported in marine invertebrate metamorphosis. The tail loss in Ciona intestinalis is regulated by an apoptosisrelated factor during larval metamorphosis [81], suggesting the necessity of apoptosis for larval metamorphosis.

Proteins associated with shell formation
Changes in shell formation are important in the development of mollusks. The oyster shell undergoes transformation during metamorphosis from an aragonite larval shell to the calcite adult shell [82]. In this study, calmodulin and calreticulin were up-regulated in the MET and eMET groups (Supplementary Table 6). Calmodulin and calreticulin were also identified in the larvae of the snail P. canaliculata during shell development [43] and in the Pacific oyster C. gigas [35]. It has been reported that calmodulin was involved in shell formation [83] and that it is associated with calcium signaling. The upregulation of both proteins suggests that they may be involved in the calcification of larval shells. Calreticulin constitutes a cellular protein "quality-control" system [84] that is important for developmental processes.
In addition, calreticulin can bind to misfolded proteins, correcting the protein folding or directing the protein toward a degradation pathway. Regarding mRNA expression levels, compared with the protein expression patterns, the mRNA expression of calmodulin exhibited similar patterns (Fig. 6), suggesting that this gene was involved in shell formation. The proteins related to calcium metabolism may be involved in calcium deposition and are likely involved in the shell formation that occurs in larvae of C. angulata.
In this study, comparison of PA and MET group, collagen alpha (GI405954419, 405,961,982) was up-regulated in the MET group. In shell matrix, collagen-like VWA domain containing proteins, fibronectin-like proteins and laminins were involved in the regulation of shell formation [85,86]. Chitin is also considered the base membrane in shell organic matrix [87]. In this study, comparison of MET and eMET group, chitinase-3-like protein 1 isoform X3 was up-regulated in the eMET group.

Nervous system remodeling during metamorphosis induced by epinephrine
The oyster foot senses a suitable substratum for settlement, suggesting biogenesis of new nerve connections. As a signal of metamorphosis, the appearance of an eyespot indicates the possible remodelling of the oyster larval nervous system. In the present study, 14-3-3 and SCO-spondin were identified. Both proteins were significantly up-regulated in the eMET group (Supplementary Table 6). The 14-3-3 protein is involved in diverse physiological processes, including growth and development, receptor signal regulation, tumorigenesis, molecular chaperone, activation or inhibition of the target protein activity [88][89][90], and innate immunity in Drosophila [91]. In this study, there were no significant differences among PL, PA, and MET groups, but the 14-3-3 protein was up-regulated in the eMET group (Supplementary Table 6). The overexpression of 14-3-3 in the body plays an important role in nerve cell differentiation [92]. The 14-3-3 proteins participate in neurodevelopment, and 14-3-3ε and ζ are associated with neurogenesis and neuronal progenitors of cell differentiation in the developing brain [93]. The 14-3-3 proteins also play important roles in left-right patterning in the axial direction of the embryo during amphibian embryogenesis [94]. In the present study, SCO-spondin was expressed at a high level in the PA group and was up-regulated in the eMET group. SCO-spondin is needed for neurogenesis during early brain development, Vera et al. [95] revealed that SCO-spondin is a key embryonic cerebrospinal diffusible factor regulating the balance between proliferation and differentiation of the brain neuroepithelial cells. The result also suggested that 14-3-3 and SCO-spondin changed in larvae of C. angulata induced by epinephrine.
Calmodulin and the 14-3-3 gamma protein are activators of tryptophan 5-monooxygenase and tyrosine 3monooxygenase [96,97], which are important for the biosynthesis of serotonin, noradrenalin and epinephrine, neurotransmitters that are key in neuronal activities [98,99]. Our study results suggested that the activation mechanisms of tryptophan 5-monooxygenase and tyrosine 3monooxygenase in C. angulata might be similar to those of mammals [96,97]. The up-regulation of 14-3-3 and calmodulin in the MET and eMET groups may indicate their roles in neuronal function and may indicate the enhanced development of the nervous system in C. angulata. This result is consistent with a report on the early larval development of the Pacific oyster C. gigas [63].
Epinephrine, dopamine, gamma-aminobutyric acid (GABA), and acetylcholine have been shown to be effective neuroendocrine compounds for inducing oyster metamorphosis [31]. In this study, neuronal acetylcholine receptor subunit alpha-10 (c88012_g1) and ryanodine receptor 44F (c96209_g1) were upregulated post epinephrine stimulation in the eSEN group. Neuronal acetylcholine receptor subunit alpha-10 is an agonist of binding that may affect all of the subunits and lead to opening of an ion-conducting channel across the plasma membrane. The channel is permeable to a range of divalent cations including calcium. One of the channels that mediates this release of Ca 2+ from the endoplasmic reticulum is the ryanodine receptor (RyR), which is essential for larval development in Drosophila melanogaster [100]. Functional assays of these key receptors will help to develop new strategies to control the settlement and metamorphosis of marine animals more effectively.

Proteins associated with cytoskeletal proteins in metamorphosis
The expression levels of cytoskeletal proteins such as paxillin, myosin regulatory light chain (MYL), vinculin and caltractin were all up-regulated in the MET stage, and the expression levels of caltractin and paxillin in the MET group were 3.9 and 8 times higher than those in PL group, respectively (Supplementary Table 6). These results indicated that skeleton proteins had undergone drastic changes during metamorphosis. Up-regulation of skeleton proteins is often found in the study of settlement and metamorphosis of marine invertebrates [40,41]. Studies have shown that skeleton proteins play an important role in the metamorphic development of amphibians [101] and insects [102]. Skeleton proteins can maintain cell structure and shape, and they provide support for cell movement, cell proliferation and cell growth. Differentially expressed proteins related to cytoskeleton regulation included integrin alpha L (ITGAL), actin-related protein 2/3 complex, subunit 5 (ARPC5) and gelsolin, which further reflected the drastic changes to the cytoskeleton during settlement and metamorphosis. The metamorphosis of oysters involves drastic tissue degradation and organ reassembly, so up-regulation of these structural proteins was inevitable.

Other proteins
In present study, the radial spoke head protein was identified, which is structural components of the axoneme and participate in producing flagellar patterns through regulation of the dynein arms [45]. Cilia-and flagella-associated protein was identified, the radial spoke head protein and cilia-and flagella-associated protein increased and then decreased from PL to MET group, this protein was up-regulated in eMET group (Supplementary Table 6), the result showed that it might be related to velum atrophy.
In this study, ferritin was identified, which increased and then decreased from PL to MET group, this protein was up-regulated in eMET group (Supplementary Table 6). Ferritin participate in neural development and is required multiple times during development [103]. Ferritin is required for embryonic and larval development [104], and also participate in shell formation [105]. The result suggest that ferritin is likely important for C. angulata embryo development.

Conclusions
In this study, the proteome of larval settlement and metamorphosis induced by epinephrine in the Fujian oyster C. angulata was studied by liquid chromatography-mass spectrometry (LC-MS). The DEPs were studied by bioinformatics analysis. The results identified the DEPs involved in the epinephrine signaling pathway, and the levels of these proteins may reflect the developmental stage of oyster larvae. These proteins were up-regulated at the stage of the pediveliger, when epinephrine treatment can induce direct metamorphosis in oysters. The proteins response of oyster Crassostrea angulata pediveliger to epinephrine was showed in Fig. 8. In a word, natural settlement and metamorphosis or epinephrine-induced direct metamorphosis involve many biological processes, and how these functions mediate the regulation process needs further research. The focus of future research will be closely related to signal pathways, and the study of protein levels will be more meaningful.

Animal manipulation
The larvae (60-75 μm) were obtained from the Zhangpu oyster aquaculture farm (Fujian Province, China), where the fertilization of C. angulata was performed. The larvae were transferred to the laboratory to continue culturing during the pediveliger stage (post eyespot larval stage). The larvae were fed daily with Isochrysis zhangjiangensis, Chlorella vulgaris, and Chaetoceros sp. The seawater temperature was maintained at 24°C-26°C, and the water was filtered using 0.22-μm filter membranes. To obtain a consistent growth rate of posteyespot larvae, 80-mesh screens (pore size: 180 μm) were used to screen immature larvae so that almost all of the individuals were in the same phase of development. After filtration, more than 95% of the larvae could be identified under a microscope as being in the pediveliger stage; 70 to 80% of the larvae used their foot to explore the substratum (Supplementary Fig. 5A). The larvae were retained in sampling boxes. Because the larvae have a preference for settlement substrates, the previous polyacetal film was replaced using a harder resin epoxy insulation board (Supplementary Fig. 5B). The individual density was about 100 ind/mL.

Sample acquisition
Oyster larvae were cultured in cement ponds. When they were at the stage of the pediveliger, the larvae were collected into 10-L buckets with mesh screens at the point when most of the larvae began to explore substrates, and the larvae were then divided into several sampling boxes. The density of larvae was about 100 Fig. 8 proteins response of oyster Crassostrea angulata pediveliger to epinephrine. Red font denote upregulated proteins, green font denote downregulated proteins ind/mL, and the larvae were at the pediveliger (PL) stage; most of the larvae attached to the bottom of the sampling box after 30 min. Some larvae that did not attach were removed with water from the sampling box, and then the bottom of the sampling box was gently washed with filtered seawater to remove the remaining non-attached larvae. Finally, the settled larvae that had attached to the sampling box were gently brushed down into filtered seawater with a soft brush and collected into freezing tubes; the seawater was removed, and the freezing tube was quickly stored in liquid nitrogen. The sample was labelled as post-attachment (PA). A subsequent PA sample was collected in another sampling box, and then new filtered seawater was added to continue culturing for 10-12 h. The larvae in this sample were collected when secondary shells were observed ( Supplementary  Fig. 5C). The sample was preserved in liquid nitrogen and labelled as post-metamorphosis sample (post-metamorphosis, MET). For each developmental stage, almost all of the collected individuals were in the same period as confirmed by microscopic examination.
To study the molecular mechanism of direct metamorphosis induced by epinephrine in Fujian oyster larvae without settlement, we also prepared samples of epinephrine-induced larvae (EPI). Coon [32] reported that the best induction effect on the pacific oyster Crassostrea gigas occurred when the concentration of epinephrine was between 10 − 4 and 10 − 5 M [32]. In this study, combined with our experimental experience, 5 × 10 − 4 M epinephrine for 2 h was used to induce larvae. After screening, about 70-80% of the pediveligers were observed to explore with their foot; this ensured that most of the larvae could be induced. The larvae were collected and transferred to the sampling box, and epinephrine was added. After about 10 s, almost all the larvae sank to the bottom of the box. After 2 h, most larvae had sunk to the bottom, but some continued to swim. After an additional 0.5 h, the larvae that sank to the bottom were collected as the epinephrineinduced larvae (EPI-sensitive larvae, eSEN) that were at the stage of post-attachment, while the larvae that had been exposed to epinephrine but continued swimming were collected as the epinephrine-insensitive larvae (EPI -insensitive larvae, einSEN). The eSEN larvae were cultured for about 10 h, at which time obvious secondary shells were observed (Supplementary Fig. 5D). The samples were labelled as post-metamorphosis induced by EPI (eMET). Each group of six groups was divided into three independent biological replicates, approximately 4000-5000 per stage, and three technical replicates were performed for each sample to ensure reproducibility.

Extraction of oyster larvae proteins
For each group sample, protein extraction solution (1 ml, 50 mM Tris 0.1% SDS, 10 mM HEPES pH 8.1) was added to the larvae samples (100 mg), and then 2-3 steel beads were added at 10 min blasts, 30 times/s, boiled at 100°C for 5 min, and then cooled to room temperature on ice. The samples were centrifuged at 17000×g for 30 min at 4°C, and the supernatant was collected. Then, samples were immediately suspended in 3-4 times volume of 20% TCA/acetone solution (w/v) with 20 mM DTT, and then placed overnight at − 20°C. The sample was then centrifuged (17,000×g, 15 min, 4°C), and the supernatant was removed. The sample was washed (0.5 mL, pre-cooled 100% acetone with 20 mM DTT) and placed on ice for five minutes, then centrifuged (17, 000×g, 15 min, 4°C). The sample was washed three times, and the supernatant was removed. The protein pellet was collected, vacuum freeze-dried, and stored at − 70°C.
Protein digestion for label-free proteomic analysis DTT buffer (1 M DTT in 8 M urea buffer per 50 μg protein) was added for a final concentration of 20 mM. Determination of protein concentration used the Bradford method. The sample was kept in a water bath (30°C) for 30 min, then cooled to room temperature. The protein samples (50 μg) were alkylated with IAA buffer (1 M IAA in 8 M urea buffer) and adjusted to a final concentration of 40 mM IAA, and the mixture was incubated in darkness for 30 min. DTT buffer was added, and the final concentration of DTT was adjusted to 10 mM. The samples then had 7 times the volume of NH 4 HCO 3 solution added to reduce the concentration of Urea to 1 M. Incubation with trypsin (1 g/50 g protein) (Promega, Madison, WI) was conducted for 12 h at room temperature. The digestion was stopped by adding 5% formic acid (the final concentration was about 0.1%).

Optimization of the desalination column
One hundred microliters of methanol was added to the desalination column; the sample was centrifuged at 100×g for 1 min, and the supernatant was removed. Then, 100 μL 70% acetonitrile, and 1% formic acid solution were added; the sample was centrifuged at 110×g for 1 min, and the supernatant was removed. Then, 200 μL 1% formic acid solution was added; the sample was centrifuged at 110×g for 1 min, and the supernatant was removed.

Sample loading and elution
Formic acid was added into the trypsin digestion solution, and the pH value was adjusted to pH 2-pH 3. The trypsin digestion solution was added into the desalination column; the sample was centrifuged at 110×g for 1 min, and the supernatant was removed. Then, 200 μl 1% formic acid solution was added; the sample was centrifuged at 110×g for 1 min, and supernatant was removed. The sample was transfered from the desalination column to a new centrifuge tube, and then 100 μl 50, 75, and 100% acetonitrile (containing 1% formic acid) were added to elute the peptide fragments from the column, and the solution was collected. The solution was then freeze-dried to powder.

Sample preparation and mass spectrometry analysis
Freeze-dried peptide segments were dissolved in 30-50 μl of a solution containing 0.1% formic acid and 2% acetonitrile. Then, the peptide segments sample was centrifuged at 12000×g for 10 min. The supernatant was filtered with a 0.22-μm microporous membrane, and the polypeptide solution was mixed in the sample bottle using an oscillator. A Thermo Scientific-Q Exactive™, a high-performance benchtop Quadrupole-Orbitrap™ Mass Spectrometer and a WATERS nano UPLC liquid phase system were employed. The Flex source was used for the ion source. The mobile phase was acetonitrile solution; phase A was 2% acetonitrile with 0.1% formic acid solution, and phase B was 98% acetonitrile with 0.1% formic acid solution.

Data analysis
Data analysis was performed with three replicate injections of each sample using a Q-Exactive mass spectrometer with the following parameters: duration: 120 min; detection method: positive ion detection; parent ion scanning scope: 300-1800 m/z. Each sample was analyzed three times. The program consisted of a 70,000 resolution full-scan MS scan, and the 10 most abundant peaks were selected for MS/MS using a 17,500 resolution scan. The ion selection window was a 1.6 massto-charge (m/z) ratio, and the normalized collision energy was 30 eV. The program used a 40s dynamic exclusion window to avoid repeated selection of peptides for MS/MS. MS1 resolution at M/Z 200: 70,000; MS2 resolution at M/Z 200: 17,500. The mass-to-charge ratios of polypeptide fragments were obtained by collecting 20 fragment patterns (MS2 scan, HCD) after each full scan.
Raw files were processed to peak lists. The original data were analyzed with three replicate injections of each sample using Proteome Discoverer software version 1.4, a protein identification software of ThermoFisher Scientific, and analyzed using SIEVE software, to quantify all of the detected peaks. The Crassostrea protein database was downloaded from NCBI (Crassostrea gigas has been sequenced with a genome-wide database), and then combined with the protein database predicted by the Crassostrea transcriptome CDS. The software Scaffold version 4.4 was selected to carry out the relative quantification of the proteins. The database search parameters were as follows: main search ppm: 6; MS/MS tolerance ppm: 20; missed cleavage: 2; De-Isotopic: ture; enzyme: trypsin; fixed modification: carbamidomethyl (C); variable modification: oxidation (M), acetyl (protein Nterm); LFQ: ture; LFQ min ratio count: 1; decoy database pattern: reserve; match between runs: 2 min; protein or peptide FDR: 0.01.

GO and KEGG annotation and enrichment
The annotations included GO and KEGG protein annotations through transcripts and the NCBI protein database based on the GI number. Fisher's Exact test was used to evaluate the significance of enrichment levels of proteins under each GO term. The protein sequences were annotated with the KEGG database. The significance of the enrichment level of each pathway was calculated by Fisher's Exact test to determine the significance level.

Experimental validation using RT-qPCR
The mRNA expression profiles of six differentially expressed genes were validated using SYBR Green fluorescent RT-qPCR with the Fujian oyster reference gene encoding elongation factor 1(EF1-F) (c233532_g1; gi|524,916,452|ref.|XP_005113003.1| PREDICTED: elongation factor 1-alpha [Aplysia californica]) and Ubiquitin (UBQ) for each sample. This SYBR Green qPCR assay was conducted using a QuantStudio 6 Flex Real-Time PCR System following the manufacturer's instructions, and relative gene expression levels were analyzed using the 2 −ΔΔCT method. This experiment was conducted in triplicate. Primer nucleotide sequences were designed using Primer Premier 6 software (Premier Biosoft, USA) (Supplementary Table 14). PCR amplification experiments were performed under the following conditions: 95°C for 30 s, then 40 cycles of 95°C for 5 s, 55°C for 30 s, and 72°C for 30 s. Three replicates were performed for each developmental stage to ensure reproducibility.