RNA-seq analysis reveals the effect of the metamorphic cue (juvenile oysters) on the Rapana venosa larvae

As a vital developmental event, metamorphosis controls the population dynamics of most marine invertebrates and affects the breeding of economic shellfish. Rapana venosa is an economically important species in China, but artificial aquaculture has hampered its metamorphosis process. Previous studies have found that juvenile oysters can effectively induce the metamorphosis of R. venosa, but the specific induction mechanism is not clear. Here, we investigated the mechanism underlying the response of R. venosa to juvenile oysters through the RNA-seq analysis. In this study, the gene set responses to metamorphosis cues (juvenile oysters) in R. venosa were identified, and GO and KEGG enrichment analyses were further performed on these gene sets. The results showed that the expression of the prototype of the class of immediate early genes, the transcription factor AP-1, was rapidly and significantly increased, and the molecular chaperone of NOS, HSP90, exhibited lower expression in the M12 group than in the control group. In contrast, the expression of inhibitors of apoptosis (IAPs) was significantly increased upon exposure to juvenile oysters. Additionally, the Wnt signaling pathway and MAPK signaling pathway were enriched in the trend analysis. These pathways may also play critical regulatory roles in the response to juvenile oysters. Taken together, the results show that competent larvae rapidly respond to the inducing effects of oysters via some immediate early genes, such as the transcription factor AP-1, which may further regulate downstream pathways such as the MAPK signaling pathway to cause subsequent changes, including a decrease in HSP90 and an increase in IAPs. These changes together may regulate the metamorphosis of R. venosa. This study provides further evidence that juvenile oysters are the metamorphosis cues of R. venosa, which may enhance our understanding of the metamorphosis mechanism in this marine invertebrate.

As a vital developmental event, metamorphosis controls the population dynamics of most marine invertebrates and affects the breeding of economic shellfish. Rapana venosa is an economically important species in China, but artificial aquaculture has hampered its metamorphosis process. Previous studies have found that juvenile oysters can effectively induce the metamorphosis of R. venosa, but the specific induction mechanism is not clear. Here, we investigated the mechanism underlying the response of R. venosa to juvenile oysters through the RNA-seq analysis. In this study, the gene set responses to metamorphosis cues (juvenile oysters) in R. venosa were identified, and GO and KEGG enrichment analyses were further performed on these gene sets. The results showed that the expression of the prototype of the class of immediate early genes, the transcription factor AP-1, was rapidly and significantly increased, and the molecular chaperone of NOS, HSP90, exhibited lower expression in the M12 group than in the control group. In contrast, the expression of inhibitors of apoptosis (IAPs) was significantly increased upon exposure to juvenile oysters. Additionally, the Wnt signaling pathway and MAPK signaling pathway were enriched in the trend analysis. These pathways may also play critical regulatory roles in the response to juvenile oysters. Taken together, the results show that competent larvae rapidly respond to the inducing effects of oysters via some immediate early genes, such as the transcription factor AP-1, which may further regulate downstream pathways such as the MAPK signaling pathway to cause subsequent changes, including a decrease in HSP90 and an increase in IAPs. These changes together may regulate the metamorphosis of R. venosa. This study provides further evidence that juvenile oysters are the metamorphosis cues of R. venosa, which may enhance our understanding of the metamorphosis mechanism in this marine invertebrate. KEYWORDS rapana venosa, metamorphic cue, RNA-Seq, juvenile oyster, trend analysis

Introduction
Rapana venosa is an economically important gastropod in China. Recently, with the increasing market demand, fishing of R. venosa has intensified, which has severely damaged wild resources. These situations have attracted widespread attention from domestic researchers, and the artificial culture industry of R. venosa will continue to develop in the future (Yu et al., 2020). However, R. venosa is a successful worldwide invader, including in the Black Sea, the Mediterranean Sea, the Adriatic and Aegean seas, the coasts of France and the Netherlands, the Chesapeake Bay on the Atlantic coast of the United States, and the Rıó de la Plata between Uruguay and Argentina, and it causes great damage to the wild resources of oysters and other bivalves in the above areas (Xue et al., 2018). This problem has also attracted widespread concern. Therefore, the invasion mechanism urgently needs to be elucidated.
Given the special ecological status of this species, in recent years, there have been a number of studies on the behavior, reproduction and development, environmental suitability, ecological effects, and technology of artificial breeding and culturing of R. venosa (Song et al., 2016a;Song et al., 2016b;Song et al., 2016c;Zhang et al., 2017;Xu et al., 2019;Yu, 2019;Yang et al., 2020a;Yang et al., 2020b;Zhang et al., 2020;He et al., 2021;Shi et al., 2022;Yang et al., 2022). However, we have found that the metamorphosis rate of R. venosa has become extremely low in its early development stage with the development of the artificial breeding and culturing industry, resulting in a low breeding survival rate (Yu et al., 2020). This low survival rate severely restricts the development of this industry.
Metamorphosis is critical for most marine invertebrates due to the sensitivity and vulnerability of organisms during this time, which affect the population dynamics and resource recruitment of marine invertebrates . Therefore, it is very important to reveal the regulatory mechanisms of metamorphosis for the development of the artificial culture industry and the illumination of the invasion mechanism. In previous studies, we have found that the metamorphosis rate of R. venosa increases by more than 60% when juvenile oysters (Crassostrea gigas, with a shell length less than 3 cm) are present, although the metamorphosis process cannot be completed spontaneously (Cavalcanti et al., 2020). Meanwhile, Xu et al. (2019) indicated that oyster reefs can significantly promote the recovery of R. venosa resources. This may partly explain why these mollusks are destroying oyster resources. However, the specific mechanism by which oysters induce metamorphosis of R. venosa larvae is not clear.
Considering the importance of metamorphosis, its regulatory mechanism has been widely researched in species including sponges, tubeworm, coral and mollusks, especially in bivalves and herbivorous gastropods (Cavalcanti et al., 2020). Most of these species use cues from bacteria. For example,  showed that in the marine sponge Amphimedon queenslandica, bacterial symbionts can play a critical role in animal development by providing their host with the arginine needed for larval settlement and metamorphosis. Shikuma et al., (2016) found that larvae of the tubeworm Hydroides e l e g an s m e t a m o r p h o s e i n r e s p o n se t o s u r f a c e -b o u n d Pseudoalteromonas luteoviolacea bacteria, further inducing the regulation of mitogen-activated protein kinase (MAPK) signaling.
Conspecific adults and food are also important cues of metamorphosis in some species. The metamorphosis of Mytilopsis sallei and Balanus amphitrite is stimulated by conspecific adults (Matsumura et al., 1998;He et al., 2021). Metamorphosis is usually food-stimulated in some gastropods, including the herbivorous gastropods Haliotis rufescens and Crepidula fornicata (Morse et al., 1979;Taris et al., 2010) and the carnivorous gastropod Onchidoris bilamellata, whose metamorphosis cues are juvenile barnacles (Rodriguez et al., 1993).
As R. venosa is a typical carnivorous gastropod, its metamorphosis cue is also its food, juvenile oysters (Yu et al., 2020), consistent with that in O. bilamellata. We have also performed many studies on the metamorphosis of R. venosa, including studies on its morphology, behavior, enzyme kinetics, symbiotic microbiota, metabolome, proteome, transcriptome and critical genes (Song et al., 2016a;Song et al., 2016b;Song et al., 2016c;Zhang et al., 2017;Yu et al., 2020;Yang et al., 2020a;Yang et al., 2020b;Shi et al., 2022;Yang et al., 2022). However, a comprehensive analysis of the metamorphosis regulation mechanism in R. venosa has not been conducted.
Transcriptomics has been widely applied to investigate the mechanisms underlying metamorphosis in invertebrates, revealing the critical pathways and genes that regulate metamorphosis and providing novel insights (Shikuma et al., 2016;He et al., 2021). Herein, to further investigate the mechanism by which juvenile oysters induce the metamorphosis of R. venosa larvae, we obtained the mRNA expression profiles of competent larvae exposed to juvenile oysters for 2 hours (M2) and 12 hours (M12) and before exposure to juvenile oysters (control, Con) via sequencing on an Illumina NovaSeq 6000 platform. Furthermore, we identified the differentially expressed mRNAs (differentially expressed genes, DEGs) and performed trend analysis of gene expression based on their transcripts per million (TPM) values to further screen out the genes that were consistently regulated during exposure to juvenile oysters. Then, these genes were subjected to Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses. The findings of this study will enhance understanding of the metamorphosis mechanism of R. venosa, which may not only aid in the development of artificial breeding and culturing but also reveal the invasion mechanism and support the protection of wild populations of oysters and other bivalves.

Larval culture and sample collection
Culturing of R. venosa larvae and sample collection were performed according to Yang et al. (2022), and the competent larvae (shell height > 1,200 mm) were used to conduct the assay. Three cement pools (3.5 × 5.2 × 1.5 m) were used to culture these larvae, and the assay was performed with a larval density of 0.1 ind/ mL. First, we randomly collected 100 larvae from each pool as the Con group. Then, the juvenile oysters (Crassostrea gigas) were placed evenly in pools (50 oysters on each scallop, added about 25000 scallop shells to each pool), and 100 larvae were collected at 2 and 12 hours thereafter from each pool as the oyster-exposed groups (M2 and M12, respectively). During the whole process of experiment and sampling, the three pools were closed with no water in or out, and maintain a stable experimental conditions (23 ± 1°C, 30‰ salinity, 7mg/L dissolved oxygen). A total of nine samples were obtained from the Con, M2 and M12 groups, which were washed with PBS and stored at -80°C for RNA extraction.

RNA extraction, library preparation and sequencing
Total RNA was extracted from the larval samples using TRIzol ® Reagent (Invitrogen, USA) according the manufacturer's recommendations, and 1% agarose gels were used to monitor the degradation and contamination of RNA. The integrity and quantification of RNA samples were evaluated using a 2100 Bioanalyzer (Agilent Technologies, USA), and RNA samples with integrity numbers (RINs) > 8.0 were used for construction of an mRNA library. The library was prepared using a TruSeq ™ RNA Sample Preparation Kit (Illumina, USA) and sequenced on an Illumina NovaSeq 6000 platform (Majorbio. Shanghai, China), and 125-bp paired-end reads were generated and deposited in the Sequence Read Ar chive (ncbi.nlm.nih.gov/sra) under BioProject PRJNA855327.

Read mapping and quantification of gene expression levels
The raw reads were trimmed and quality-controlled via SeqPrep and Sickle with default parameters, and the clean reads were separately aligned to the reference genome of R. venosa (unpublished data) with orientation mode using HISAT2 software (Kim et al., 2015). The mapped reads of each sample were assembled by StringTie via a reference-based approach (Pertea et al., 2015), and then the paired-end clean reads were aligned to the reference genome of R. venosa (unpublished data) using TopHat v2.0.12. HTSeq v0.6.1 was used to count the reads mapped to each gene. Known genes and novel transcripts were identified by Reference Annotation Based Transcript (RABT) assembly method, the Cufflinks v2.1.1. Subsequently, the TPM of each gene was calculated based on the length of the gene and the number of reads mapped to the gene.

Differential expression analysis and trend analysis of genes
DEG analysis was conducted among the three groups, and the TPM algorithm was used to normalize the expression level of each transcript. RSEM was used to quantify gene abundances. Briefly, DEG analysis was performed using the DESeq2 R package with the significance thresholds set to a P-adjusted value < 0.05 and a | log 2 FC| ≥ 1. Additionally, after filtering out the genes with TPM = 0 and normalizing the TPM by log2 transformation, the trend analysis of R. venosa transcriptome data was performed on the Majorbio Cloud Platform by Short Time-series Expression Miner (STEM) software.

Functional enrichment analysis
To further explore the biological functions of the DEGs and these genes with significant trends, GO and KEGG enrichment analyses were performed. GO enrichment analysis was implemented using the GOseq package in R, and KOBAS was used to analyze the statistical enrichment of genes in KEGG pathways. GO terms and KEGG pathways with the significance enrichment threshold set to a Padjusted value < 0.05 corrected by the Bonferroni method.

Verification by quantitative real-time PCR (qRT−PCR)
To validate the accuracy of the RNA-Seq profiling results, critical genes were selected for qRT−PCR analysis. The cDNA for qRT−PCR was synthesized using the Prime Script ™ RT Reagent Kit with gDNA Eraser (TaKaRa, Japan). The primers used in the mRNA qRT−PCR assay were designed using Primer 5, and 60S ribosomal protein L28 (RL28) was selected as a housekeeping gene to normalize the data . A SYBR PrimeScript RT−PCR Kit II (TaKaRa, Japan) was used to quantify the expression levels. The relative expression levels of genes were estimated by the 2 −DDCT method. All data are presented as the means ± SEs (N = 3). Statistical significance was analyzed using SPSS v.19, with a P value < 0.05 considered to indicate significance.

Overview of RNA-seq and quantitative RT-PCR validation
Nine libraries of R. venosa larvae were constructed from the Con group and oyster-exposed groups (M2 and M12), from which 63.84 Gb of clean reads were generated. Each cDNA library was greater than 6.09 Gb in size (Q30 > 93.36%), and the rate of successful mapping to the reference genome ranged from 75.32 to 76.71% (Table 1). A total of 37,493 genes were obtained, including 24,662 known genes and 12,831 novel genes. Gene function annotations showed that 19,629 genes (79.59% of known genes) had significant matches in the COG, GO, KEGG, KOG, Pfam, SwissProt, eggNOG, or NR databases (Table 2). Nine random genes were then chosen to verify the accuracy of the transcriptome sequence data by qRT-PCR. The results were shown in Figure 1. Altogether, the transcriptome data were consistent with the qRT-PCR results, which further reinforced our findings.

Sample clustering and heatmap analysis
Hierarchical clustering analysis (HCA) showed that the groups exposed to oysters for 12 hours (M12) were clustered as one branch, while the group exposed to oysters for 2 hours (M2) was clustered with the Con group ( Figure 2A). Additionally, principal component analysis (PCA) showed that all three groups were significantly separated from each other ( Figure 2B). These results may mean that exposure to oysters has a significant effect on the mRNA profile of R. venosa larvae and that the effect increases with time.

Differential expression analysis and functional enrichment analysis
Differential expression analysis of genes among the three groups was also performed based on the TPM value. The DEGs in the Con vs. M2, Con vs. M12 and M2 vs. M12 comparisons were regarded as those induced by juvenile oysters but participating in different time periods of the induced-regulation process. In total, there were 189 DEGs in Con vs. M2, 597 DEGs in Con vs. M12 and 744 DEGs in M2 vs. M12 ( Figure 3A). Additionally, a Venn diagram showed that 42 DEGs were shared between the Con vs. M2 and Con vs. M12 comparisons, 55 DEGs were shared between the Con vs. M2 and M2 vs. M12 comparisons, and 292 DEGs were shared between the Con vs. M12 and M2 vs. M12 comparisons ( Figure 3B). In contrast, 92 DEGs were unique to Con vs. M2 comparisons, 263 DEGs were unique to Con vs. M12 comparisons, and 397 DEGs were unique to M2 vs. M12 comparisons.
The 189 DEGs in the Con vs. M2 comparison were regarded as the initial-response genes (I-DEGs) to oyster presence, which were significantly enriched with 2 GO terms: DNA-binding transcription factor activity and transcription regular activity ( Figure 4A). There were no significantly enriched pathways (P adjust < 0.05). We further analyzed the expression of these genes enriched with the 2 GO terms. As shown in Figure 4B, 12 genes were upregulated in the M2 group compared with the Con group, including some nuclear receptor genes (nuclear receptor subfamily 1 and nuclear hormone receptor E75), CREBRF, the transcription factor Maf/MafK/AP-1, dendritic arbor reduction protein 1, cAMP-responsive element-binding protein, CCAAT/enhancer-binding protein beta and transforming protein v-Fos/v-Fox. Five other genes were downregulated in the M2 group, including the transcription factor MafAa, forkhead box protein B1, vitamin D3 receptor, deformed epidermal autoregulatory factor 1 and nuclear factor interleukin-3-regulated protein.
The 597 DEGs in Con vs. M12 were regarded as the final-response genes (F-DEGs) to oyster presence, which were significantly enriched in 3 pathways, including ribosome biogenesis in eukaryotes, antigen processing and presentation, and RNA transport (P adjusted < 0.05) ( Figure 5A). There were 36 DEGs in the 3 pathways, and all of them were downregulated in the M12 group compared with the Con group, which suggested that these 3 pathways were suppressed in the presence of oysters or were inhibited by oyster exposure. There were 36 DEGs in the 3 pathways, including many heat shock  proteins (HSPs), such as HSP70 proteins (heat shock protein 68, heat shock 70 kDa protein A1, heat shock cognate 71 kDa protein and so on) and HSP90 proteins (97 kDa heat shock protein and heat shock protein HSP 90-beta) ( Figure 5C). These DEGs were significantly enriched with 52 GO terms (top 20 were shown), and the smallsubunit processome term was the most enriched ( Figure 5B). In addition to the initial-response and final-response genes that were altered in response to oysters, 397 DEGs existed only in the M2 vs. M12 comparison. These genes may respond to oyster presence slowly and ultimately not significantly change but still have an important regulatory effect; therefore, they should not be ignored. We regarded them as late-response genes (L-DEGs) to oyster exposure. The 397 DEGs were significantly enriched with 3 GO terms, including rRNA metabolic process, ncRNA metabolic process and rRNA processing. Furthermore, we analyzed the expression of 14 genes enriched with 3 GO terms ( Figure 6A). All 14 genes were downregulated in the M12 group compared with the M2 group, which may mean that the metabolic processes related to rRNA and ncRNA were slowed down in the later stage in the presence of oysters ( Figure 6B).  Validation of the 9 random DEGs involved in different pairwise comparisons by qRT-PCR. RL28 was selected to normalize the gene expression levels. The data are shown as means ( ± SE) of three replicates. Different superscripted letters indicate significant differences (P < 0.05).

Analysis of gene expression trends and functional enrichment analysis
A total of 22,662 genes were used to perform trend analysis, and 8 model profiles were used to summarize the expression patterns. The genes in the same profile presented similar expression patterns. As shown in Figure 7A, 3 profiles labeled with different colors were statistically significant, including profile 2 (3563 genes), profile 4 (1308 genes), and profile 1 (751 genes). These genes with different expression trends may have different functions in the metamorphosis induced by oysters. To further understand the functions of these genes in the 3 profiles, GO and KEGG enrichment analyses were performed. The top 20 significantly enriched GO terms and pathways are illustrated. For the genes in profile 2, 10 of the significantly enriched pathways were in the "human disease" category, 5 were in the "organismal systems" category, and only one was in the "cellular process" category. Four were in the "Environmental Information Processing" category, including the Ras signaling pathway (62 genes with adjusted P < 0.05), the Wnt signaling pathway (46 genes with adjusted P < 0.05), the MAPK signaling pathway (62 genes with The functional enrichment and expression analysis of initial-response genes to oyster induction. (A) GO enrichment results of 92 DEGs in Con vs. M2 comparison. (P adjust < 0.05). (B) Heat map of the expression analysis of 17 DEGs in GO:0003700 (DNA-binding transcription factor activity) and GO:0140110 (transcription regulator activity). Yang et al. 10.3389/fmars.2023.1122668 Frontiers in Marine Science frontiersin.org adjusted P < 0.05), and the Rap1 signaling pathway (64 genes with adjusted P < 0.05) ( Figure 7B). In addition, the genes in profile 2 were significantly enriched in 10 GO terms belonging to the "biological process" category, 1 belonging to the "cellular component" category, and 9 belonging to the "molecular function" category, which also contained the Wnt signaling pathway (22 genes with adjusted P < 0.05) ( Figure 7C). For the genes in profile 4, only one pathway, autophagy -animal (18 genes with adjusted P < 0.05), was significantly enriched, and no significantly enriched GO term was identified. For the genes in profile 1, there were also no significantly enriched GO terms or pathways. Furthermore, we analyzed the expression of the genes in the Wnt signaling pathway and animal autophagy.

Discussion
The food or prey of some animals not only provides them with the necessary nutrients for growth and development but also plays an important role in their development at a critical time. One of the most dramatic examples is when the food or prey induces animal metamorphosis (Morse et al., 1979;Rodriguez et al., 1993;Taris et al., 2010). Although representative species from diverse taxa undergo food-induced metamorphosis, little is known about the processes that occur within the animals as they detect the cues produced by their food. To further understand the mechanism by which juvenile oysters induce the metamorphosis of R. venosa larvae, we performed RNA-seq and analyzed changes in gene expression at different time points upon exposure to metamorphosis cues, the juvenile oysters. Then, we conducted functional enrichment for these genes and identified the critical signaling pathways in the process. The present results may provide new insights into the regulatory mechanisms of metamorphosis in marine invertebrates.
In the present study, GO enrichment analysis revealed that 189 DEGs that responded early to the presence of juvenile oysters were significantly enriched with two GO terms: DNA-binding transcription factor activity and transcription regulator. These genes included some critical transcription factor genes, such as Maf, AP-1, MafK, MafAa, and transforming protein v-Fos/v-Fox, and most of the genes were upregulated. A previous study has also indicated that these genes are Jun or Fos genes (transcription factor AP-1 gene family), which are prototypical immediate early genes. Immediate early genes are characterized by a rapid and transient activation of transcription in response to changes of environmental conditions, such as growth factors, cytokines, tumor promoters, carcinogens, and expression of certain oncogenes (Angel and Hess, 2010). The transcription factor AP-1 gene family participates in many pathways, such as the MAPK signaling pathway, Wnt signaling pathway and apoptosis, which have been proven to be important in the metamorphoses of marine invertebrate animals (Angel and Hess, 2010). Therefore, we speculate that the transcription factor AP-1 may be a critical initiator in competent larvae of R. venosa in response to the metamorphosis cues of juvenile oysters that regulates other subsequent processes.
A previous study has indicated that the most critical members of the class of protein kinases regulating the activity of AP-1 in response to extracellular stimuli are MAPKs (Angel and Hess, 2010). In addition, researchers have shown that the MAPK signaling pathway is induced by the metamorphosis cue in tubeworm Hydroides elegans and regulates the downstream process of metamorphosis (Shikuma et al., 2016;Wang and Qian, 2010), which also occurs in other marine invertebrate animals, such as corals, ascidians and barnacles (Chambon et al., 2007;He et al., 2012;Siboni et al., 2014). Interestingly, we also found that the MAPK signaling pathway was significantly enriched in the gene set of profile 2 in the trend analysis, which may mean that in competent larvae of R. venosa, MAPK signaling participates in the response to metamorphosis cues. Chambon et al. (2007) found that the MAPK signaling pathway regulates gene networks that stimulate metamorphosis and apoptosis in tail tissues of ascidian tadpoles, and apoptosis occurs during metamorphosis in almost all invertebrates, including R. venosa (Song et al., 2016a). In the present study, we found that the expression of IAPs was significantly increased upon exposure to juvenile oysters. IAP expression is also upregulated after metamorphosis in juveniles of R. venosa (Shi, 2022). IAPs are important regulators that participate in the sophisticated molecular machineries in cells that guard against inappropriate or premature apoptosis to counter various stimuli capable of triggering death . IAPs are also highly upregulated in response to various environmental stresses, including heat, air exposure, and low-salinity stress, in some marine invertebrate animals . The presence of oysters can also be considered an environmental factor for competent R. venosa larvae, and the upregulation of IAPs in the present results indicates that apoptosis is tightly regulated during the metamorphosis of R. venosa. The observed upregulation also suggests that apoptosis is induced by the presence of oysters during the metamorphosis of R. venosa. Furthermore, we found that the transcription factor AP-1 affected the expression of some proapoptotic genes, including p53, Fas, Fas-L, Bim and HRK, which supports the findings of a previous study showing that AP-1 target genes are involved in the regulation of apoptosis (Angel and Hess, 2010).
Past studies have revealed that, as a molecular chaperone of nitric oxide synthase (NOS), HSP90 is also a negative regulator (repressor) of metamorphosis in a diverse range of marine invertebrates (Bishop and Brandhorst, 2001;Bishop et al., 2001;Hens et al., 2006;Taris et al., 2009). However, Ueda and Degnan (2014) had a different discovery: that inhibition of NOS reduces rates of metamorphosis, which means that NO facilitates, rather than represses, induction of metamorphosis in Haliotis asinina. Additionally, Ueda and Degnan (2014) tested the hypothesis of NOS-HSP90 interaction at the initiation of settlement and metamorphosis because they found that the spatial and temporal expression trends of NOS and HSP90 were not consistent. In the present study, we found that HSP90 genes were significantly downregulated after exposure to oysters for 12 hours ( Figure 4C), which might be considered to be a signal of metamorphosis induced by the oysters, but the specific function of HSP90 needs further investigation.
In addition to HSP90 genes, several HSP70 genes were also significantly downregulated in oysters ( Figure 4C). Similar trends of downregulation were observed in a study on gene expression during the metamorphosis of R. venosa . HSP70 first received attention for its role in the response to heat shock (Welch, 1993) and has also been found to be critical for the folding and assembly of other cellular proteins (Ritossa, 1962). HSP70 genes also participate in the regulation of the MAPK signaling pathway (Patel, 2009), which has been widely reported to participate in the regulation of metamorphosis in several marine invertebrate animals mentioned above. Porto et al. (2018) reported that HSP70 can act as a signaling molecule, modulating MAPK downstream signaling during memory consolidation in the hippocampus. Patel (2009) also indicated that HSP70 plays a role in the induction of p38MAPK in animals to respond differently to high osmotic stress. Therefore, we speculate that HSP70 may also participate in the MAPK signaling pathway and be related to the metamorphosis of R. venosa; however, the specific effect of HSP70 on metamorphosis needs further investigation, especially since no previous studies have reported the relationship between HSP70 and metamorphosis in marine invertebrates.
The Wnt signaling pathway, which has been widely studied, regulates cell proliferation, migration, differentiation and survival during animal development as well as the maintenance of homeostasis and regeneration (Liu et al., 2008;De Robertis, 2010;Clevers and Nusse, 2012). In the trend analysis, we found that the Wnt signaling pathway was significantly enriched in profile 2, a gene set that was decreased at the early stage during exposure to juvenile oysters. This finding may suggest that the Wnt signaling pathway participates in the metamorphosis of R. venosa. Chandramouli et al. (2013) also found Wnt signaling pathway enrichment during the larval metamorphosis of the polychaete worm Pseudopolydora vexillosa in transcriptome studies. Furthermore, Yue et al. (2012) have reported that Wnt signaling pathways may be involved in the metamorphosis of the bryozoan Bugula neritina. However, the specific function of the Wnt signaling pathway in the metamorphosis of marine invertebrates needs further investigation.

Conclusion
In this study, the gene set responses to metamorphosis cues (juvenile oysters) in R. venosa (I-DEGs and F-DEGs) were identified, and GO and KEGG enrichment analyses were further performed on these gene sets. We found that DNA-binding transcription factor activity and transcription regulator activity were enriched for the I-DEGs, such as the transcription factor AP-1, which is a prototypical immediate early gene that rapidly and transiently activates transcription in response to changes in environmental conditions. The molecular chaperone of NOS, HSP90, exhibited lower expression in the M12 group than in the Con group, while the expression of IAPs was significantly increased upon exposure to juvenile oysters. Additionally, the Wnt signaling pathway and MAPK signaling pathway were enriched in the trend analysis. These pathways may also play critical regulatory roles in the response to juvenile oysters. Taken together, the results show that competent larvae rapidly respond to the presence of oysters by expressing some immediate early genes, such as the transcription factor AP-1. These genes may further regulate downstream pathways such as the MAPK signaling pathway and result in subsequent changes, including a decrease in HSP90 and an increase in IAPs. These changes together may regulate the metamorphosis of the larvae.

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics statement
The studies involving animals were reviewed and approved by the Animal Welfare Committee of the Institute of Oceanology, Chinese Academy of Sciences.   Science and Technology (no. LMEESCTSP-2018). The funders had no role in the study design, data collection and analysis, decision to publish or preparation of the manuscript.

Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.