Understanding microRNA Regulation Involved in the Metamorphosis of the Veined Rapa Whelk (Rapana venosa)

The veined rapa whelk (Rapana venosa) is widely consumed in China. Nevertheless, it preys on oceanic bivalves, thereby reducing this resource worldwide. Its larval metamorphosis comprises a transition from pelagic to benthic form, which involves considerable physiological and structural changes and has vital roles in its natural populations and commercial breeding. Thus, understanding the endogenous microRNAs (miRNAs) that drive metamorphosis is of great interest. This is the first study to use high-throughput sequencing to examine the alterations in miRNA expression that occur during metamorphosis in a marine gastropod. A total of 195 differentially expressed miRNAs were obtained. Sixty-five of these were expressed during the transition from precompetent to competent larvae. Thirty-three of these were upregulated and the others were downregulated. Another 123 miRNAs were expressed during the transition from competent to postlarvae. Ninety-six of these were upregulated and the remaining 27 were downregulated. The expression of miR-276-y, miR-100-x, miR-183-x, and miR-263-x showed a >100-fold change during development, while the miR-242-x and novel-m0052-3p expression levels changed over 3000-fold. Putative target gene coexpression, gene ontology, and pathway analyses suggest that these miRNAs play important parts in cell proliferation, migration, apoptosis, metabolic regulation, and energy absorption. Twenty miRNAs and their target genes involved in ingestion, digestion, cytoskeleton, cell adhesion, and apoptosis were identified. Nine of them were analyzed with real-time polymerase chain reaction (PCR), which showed an inverse correlation between the miRNAs and their relative expression levels. Our data elucidate the role of miRNAs in R. venosa metamorphic transition and serve as a solid basis for further investigations into regulatory mechanisms of gastropod metamorphosis.

undertaken by many enterprises since 1992 owing to its economic importance (Yuan 1992). However, large-scale aquaculture of this species has been hampered by difficulties with larval culture during metamorphosis. Moreover, wild veined rapa whelk resources have been declining in China owing to increasing fishing activity. This whelk is considered an invasive species beyond the western Pacific Ocean as a result of unintended world-wide transportation, and it heavily threatens the biomass of local bivalves (Yuan 1992). It was first recorded as an invasive species in the Black Sea during the 1940s (Drapkin 1963). Since then, primarily owing to unintended global transport, R. venosa has become extremely pervasive and has extended its range to Quiberon Bay, France (Mann et al. 2004); Chesapeake Bay (Harding and Mann 1999); Rio de la Plata between Uruguay and Argentina (Pastorino et al. 2000); and The Netherlands' coastal waters (Nieweg et al. 2005). Its prevalence heavily disrupts native trophic structure and damages endemic bivalve resources. Metamorphosis may control population dynamics; therefore, understanding its mechanism may be useful for aquaculture, resource restoration, and preventing biological invasion of R. venosa.
There are few published reports on R. venosa metamorphosis. A previous study described morphological changes that R. venosa undergoes in the metamorphosis process; the study also indicated that during this process, the diet of R. venosa shifts from phytophagous to carnivorous (Pan et al. 2013). Metamorphosis inducers for this species were also investigated, and it was found that acetylcholine chloride and calcium chloride (CaCl 2 ) were effective and had low toxicity (Yang et al. 2015), suggesting that these compounds may be promising in artificial aquaculture. A comprehensive transcriptome database of R. venosa has been constructed from precompetent, competent, and postlarvae (Song et al. 2016c), and it forms a baseline for future studies on gene/protein activity associated with metamorphosis. Transcriptomic and proteomic analyses of R. venosa metamorphosis identified differentially expressed genes/proteins. This finding indicates that there are multiple processes involved in its metamorphic transition, especially ingestion and digestion, cytoskeleton and cell adhesion, stress response and immunity, and tissue-specific development (Song et al. 2016b,d). The metabolic profiles of competent and postlarval stages of R. venosa were examined by gas chromatography-mass spectrometry (GC2MS). The analysis detected 53 metabolites whose concentrations differed before and after metamorphosis. They are indicative of the changes in energy metabolism and cell signaling that occur during metamorphosis (Song et al. 2016a). Nevertheless, since microRNAs (miRNAs) participate in RNA silencing and post-transcriptional gene expression regulation (Ambros 2004;Bartel 2004), miRNA data are required to provide further concrete support for conclusions drawn from the transcriptome/proteome data.
The miRNAs are short ($22 nt) and noncoding and have been implicated in cell differentiation, proliferation, migration, and apoptosis. They downregulate the expression of target genes by binding to their 39 untranslated regions (UTRs) (Sun et al. 2017). Metamorphosis is essential for developmental and evolutionary success. Nevertheless, the mechanism by which miRNAs regulate this process remains to be determined. Over the past 10 yr, several studies have been conducted on this aspect in insects, fish, and amphibians. A mutation was found that eliminates let-7 and miR-125 and leads to widespread defects during the metamorphosis of the insects Blattella germanica and Drosophila melanogaster (Caygill and Johnston 2008;Mercedes et al. 2012;Sempere et al. 2002). When anti-miR-100 depleted miR-100 in the last instar of the hemimetabolan insect B. germanica, the adult wings were slightly smaller than those of the wild type. They also presented with partial fusion of the cross-veins in the anterior remigium and abnormal bifurcations of those in the posterior remigium. The depletion of let-7 elicited the same adult wing vein pattern malformations (Rubio and Belles 2013). The MiR-2 family regulates B. germanica metamorphosis by controlling the juvenile hormone signaling pathway (Lozano et al. 2015). In the holometabolan D. melanogaster, miR-9a-mutants showed wing margin defects and a few ectopic sensory organs (Li et al. 2006), and in a later study, miR-9 was found to prevent apoptosis during wing development (Bejarano et al. 2010). In the Japanese flounder (Paralichthys olivaceus), the expression patterns of 197 miRNAs during metamorphosis were analyzed, and a later study showed that the decrease in miR-17 upregulated Cdc42 during metamorphosis (Zhang et al. 2016). In the amphibian Xenopus tropicalis and in fish, expression-profiling miRNAs at metamorphosis were identified. The miR-133 played an important part in skeletal muscle development during metamorphosis (Zhang et al. 2016). Molluscs include the largest marine phylum, and comprise $23% of the total marine organisms. Metamorphosis is the most important developmental event in the molluscan life cycle; however, the characterization and roles of miRNAs in molluscan metamorphosis have not been determined to date.
The purpose of the present study was to elucidate the endogenous miRNAs that drive metamorphosis in the veined rapa whelk R. venosa. By sequencing on the Illumina HiSeq 2500 platform, we compared the global expression profiles of small RNAs in precompetent larvae (pre-CL), competent larvae (CL), and postlarvae (PL). In previous studies, we investigated the mRNA global expression profile of whelk metamorphosis (Song et al. 2016d); therefore, in the present study we performed a differentially expressed miRNA-mRNA correlation analysis to elucidate miRNA regulation in whelk metamorphosis. These findings will provide new insights into gastropod metamorphosis and facilitate investigation of miRNA function in a biphasic life cycle in the future.

Sampling
Parent R. venosa were collected from their naturally growing areas in Laizhou Bay, China (37°1194.78$N, 119°4193.75$E). Parent culture, spawning, larval incubation, and rearing were performed at the Blue Ocean (Laizhou, Shandong, China) according to previously published methods (Pan et al. 2013). Planktonic larvae were cultured in 2.5 · 4 · 1.2 m cement pools at 23.5-25.8°and a density of 0.3/ml. Isochrysis galbana, Chlorella vulgaris, and Platymonas subcordiformis were pooled and provided as a daily diet (2 · 10 5 cells/ml, two times) to the pelagic larvae. Samples from the precompetent larval (three spiral-whorls) stage, the competent larval (four spiral-whorls) stage, and the postlarval (juvenile) stage were collected as three biological replicates, each of which consisted of 40-100 individuals. The samples were inspected under a microscope to ensure that .95% of the individuals were developmentally synchronized. Each replicate was then rinsed with doubledistilled H 2 O and flash-frozen in liquid nitrogen until use.

Library construction and sequencing
Total RNA was extracted from an individual intestine using the RNeasyMini Kit (Qiagen, Germantown, MD) according to the manufacturer's instructions. The quality and concentration of RNA were measured using a NanoDrop 1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA). RNA molecules in the size range of 18-30 nt were enriched by polyacrylamide gel electrophoresis (PAGE). The 39 adapters were added and the 36-44 nt RNAs were enriched. The 59 adapters were then ligated to the RNAs as well. The ligation products were reversetranscribed by polymerase chain reaction (PCR) amplification. The 140-160-bp PCR products were enriched to generate a cDNA library and sequenced using Illumina HiSeq 2500 (Gene Denovo Biotechnology, Guangzhou, China).

Sequence data analysis
Raw data were obtained from base calling on the original image. They were cleaned by removing reads containing poly-N or poly-A/T/C/G and 59 adapter contaminants, and those missing 39 adapters or insert tags. Low-quality reads were also eliminated. Other RNAs (tags originating from protein-coding genes, repeat sequences, rRNA, tRNA, snRNA, and snoRNA) were deleted after blasting against the Repeat-Masker (www.repeatmasker.org/), the GenBank database (http://blast. ncbi.nlm.nih.gov), and the Rfam database (http://sanger.ac.uk/software/ Rfam). No miRNA information for rapa whelk is included in miRBase v. 21.0, so the remaining clean reads were aligned to it with #2 mismatches to seek all known precursor/mature miRNAs. The miRNAs with the highest expression for each mature miRNA family were selected as temporary miRNA references. Clean data were aligned to them and their expression levels were calculated by summing the read counts aligned with the temporary miRNA database with #2 mismatches. The precursors of the identified miRNAs were predicted. Molecules without a hairpin structure were identified as pseudo-miRNAs. The potentially novel miRNAs were detected by using MIREAP (http://sourceforge.net/projects/mireap/) with stemloop structure prediction (Chen et al. 2013).
By pairwise comparison of the miRNA expressions among precompetent, competent and postlarval samples, the differentially expressed miRNAs were identified and estimated. This procedure conforms to the BGI standard as follows: (1) the expression of the miRNA in each sample was normalized to determine the TPM (transcripts per million). The formula TPM = actual miRNA count/total count of clean reads · 10 6 was normalized.
(2) The final TPM in each larval stage was calculated by averaging three biological replicates. (3) Fold-change and P-values were calculated based on the normalized expression as described previously (Sun et al. 2017).

Quantitative real-time PCR
The total RNA isolated as described above was used in real-time PCR. Briefly, 1 mg of total RNA was reverse-transcribed into cDNA using the One Step PrimeScript miRNA cDNA Synthesis Kit (TaKaRa Bio, Kusatsu, Shiga, Japan), following the manufacturer's directions. For the mRNA quantitative PCR (qPCR) assay, the primers (Supplemental Material, Table S1) were designed using Primer3 (http://primer3. sourceforge.net/releases.php). The reaction proceeded as follows: 95°f or 3 min; 95°for 15 sec, 60°for 25 sec, 72°for 15 sec for 40 cycles, and a final 20-min extension step at 72°. RL28 was chosen as a control gene for internal standardization (Song et al. 2017). The mRNA expression levels were quantified with the SYBR PrimeScript RT-PCR Kit II (TaKaRa Bio) following the manufacturer's instructions on an Eppendorf Mastercycler ep realplex platform (Eppendorf, Hamburg, Germany). For the miRNA qPCR assay, the reaction proceeded for 60 min at 37°and 5 sec at 85°. The cDNA was amplified using realtime PCR and platinum SYRB Green qPCR SuperMix-UDG (Invitrogen, Carlsbad, CA) with miRNA-specific forward and reverse primers (Table S1). The 5.8 S rRNA was used as an internal reference gene to normalize the data (Zheng et al. 2016). The amplification products were detected by melting curve and gel electrophoresis to ensure primer efficiency and PCR specificity. The relative expression levels of mRNA and miRNA were estimated by the 2 2DDCT method. All data were presented as means 6 SE (N = 3). Statistical significance was analyzed using SPSS v. 18, with P , 0.05 being considered significant.  (Song et al. 2016d) and the differentially expressed miRNAs were integrated to analyze the key miRNA-target pairs. Only the inversely correlated miRNA-target pairs with MFE # 218 were screened.

Data availability
Raw sequencing data were submitted to the GEO (Gene Expression Omnibus) database with accession No. GSE102631. Supplemental materials include the details of primers for qPCR assays (Table S1), the differentially expressed miRNAs during metamorphosis development (Table S2), and miRNA-target pairs of differentially expressed miRNAs (Table S3).

MicroRNA library construction
To identify the miRNAs differentially expressed during metamorphosis in rapa whelk, nine small RNA libraries (precompetent larvae: Pre-CL 1, Pre-CL 2, and Pre-CL 3; competent larvae: CL 1, CL 2, and CL 3; postlarvae: PL 1, PL 2, and PL 3) were constructed using Illumina sequencing. A dataset consisting of $85,000,000 reads (ranging from 8,235,539 to 10,159,366 among the samples) was obtained after trimming the adapter sequences (

Selection of miRNA-target pairs and qPCR validation
The aforementioned mRNA expression profiling data from the same metamorphosis sampling stages (Song et al. 2016c) were used to perform association analyses along with the current miRNA profiling. The differentially expressed miRNAs detected in the present study were used to select the miRNA-target pairs expressed in metamorphosis (Table S3). The miRNAs negatively regulate the expression of their target mRNAs either by translation inhibition or by mRNA degradation. In previous studies (Song et al. 2016b,c), the mRNAs/proteins involved in "ingestion and digestion," "cytoskeleton and cell adhesion," and "apoptosis" were thought to have important roles in driving metamorphosis. We identified 20 key miRNA-target pairs potentially implicated in these aspects of whelk metamorphosis (Table 3). For example, in "ingestion and digestion," we found that let-7-y potentially regulates the SARP-19 precursor, conotoxin Cl14.12, and the exoglucanase XynX genes. Therefore, a single miRNA may regulate multiple target genes during metamorphosis. The miR-1175-x targets the cysteine-rich secretory protein gene, and miR-2001-x targets endo-1,4-b-xylanase. The gene miR-71-x regulates the b-1,4-xylanase and membrane metalloendopeptidase-like 1 genes. In "cytoskeleton and cell adhesion," tektin-3, dynein heavy chain 8 (axonemal), and dynein intermediate chain 2 (ciliary), all of which are the main structures of velum cilia, were regulated by miR-5106-y, miR-87-y, and miR-315-x, respectively. Novel-m0020-5p determines the expression of n apoptosis 2 inhibitor genes, and miR-276-y regulates the caspase-3 gene. These two mRNAs are both important for programmed cell death during metamorphosis.
Real-time PCR analysis was performed on nine key targeted miRNA-mRNA pairs (Figure 4) to validate and identify the metamorphosisrelated miRNAs in rapa whelk. The results showed that the miRNAs were consistent with the overall trend in high-throughput sequencing. For each of the nine pairs, there was an inverse correlation between the expression levels of miRNA and miRNA. For example, the let-7-y miRNA decreased during the transition from precompetent to postlarval, whereas its target mRNAs (c112229_g1 SARP-19 precursor, c124801_g1 Conotoxin Cl14.12, and c150903_g1 Exoglucanase XynX) significantly increased.

DISCUSSION
In this study, we investigated the miRNA expression profile during the metamorphosis of the rapa whelk using high-throughput sequencing. In total, nine libraries were constructed and 85,000,000 reads were obtained. These results will augment information on the small RNA genome of rapa whelk and provide a basis for miRNA regulation during metamorphosis. A total of 195 miRNAs was significantly differentially expressed among three larval stages (precompetent, competent, and postlarval stages) and targeted thousands of genes. This result was expected since metamorphosis is the most complicated developmental event of the life cycle and entails considerable structural, physiological, and behavioral changes.
To obtain an insight into the possible functions of the differentially expressed miRNAs involved in metamorphosis, we performed GO and KEGG pathway enrichment analyses on their predicted target genes. The significantly enriched GO terms "biological adhesion," "cell aggregation," "cellular component organisation or biogenesis," "localisation," "developmental process," "signalling," "immune system," and "response to stimulus" are primarily associated with development, gene expression, immunity, and the cell cycle. Metamorphosis is a complex process and includes tissue remodeling, cell migration, differentiation, proliferation, and others (Jackson et al. 2005). A total of eight significantly enriched pathways were observed. The "TNF signalling" pathways were enriched because they trigger apoptosis, and old organs such as the velum are degenerated by programmed cell death. "SNARE interactions in vesicular transport" were enriched because the nervous system mediates metamorphosis in many marine invertebrates (Voronezhskaya 2004;Gifondorwa and Leise 2006) and SNARE participates in vesicle docking, priming, fusion, and the synchronization of neurotransmitter release into the synaptic cleft during neurosecretion (Shi et al. 2011). SNAREs also play a crucial part in the autophagy required for velum degradation and reabsorption during molluscan metamorphosis. "Nicotinate and nicotinamide metabolism," "Ubiquitin mediated proteolysis," "Phosphonate and phosphinate metabolism," "Pyrimidine metabolism," and "Sulphur relay system" were also affected because of the tissue remodeling and energy redistribution that occur during metamorphosis. Table 2 lists 39 differentially expressed miRNAs with striking changes. The miR-242-x expression level steadily decreases as precompetent larvae develop into postlarvae. Its expression level decreased by .3000 fold after metamorphosis. Therefore, the miR-242 family may have important roles Figure 3 Top 20 pathway enrichment of the predicted target genes of differentially expressed miRNAs. The enrichment factor is the ratio of the number of target genes of differentially expressed miRNAs annotated in this term to the number of all genes annotated in the same term.
in rapa whelk metamorphosis. In Caenorhabditis elegans, miR-242 and miR-793 target the Argonaute protein ALG-1, which controls the RNA interference process involved in developmental timing (Grishok et al. 2001). The expression level of novel-m0052-3p in precompetent larvae remained at an average TPM of 8.17. On the other hand, the average TPM decreased to 0.01 in the competent larvae and rose sharply to 36.61 in the postlarvae. The function of this miRNA remains as yet unknown and, to our knowledge, no relevant studies on it have been published to date. The TPM level of miR-100-x in the precompetent larvae was 59.29. It decreased to 23.12 in the competent larvae but sharply increased to 3974.4 after metamorphosis. There was a similar trend for miR-125-x. Both miR-100 and miR-125 are believed to participate in cell migration. Low levels of miR-100 and miR-125 may promote hepatocellular carcinoma metastasis (Rubio and Belles 2013). Rapa whelk metamorphosis is accompanied by high levels of cell death, proliferation, and tissue remodeling, thus involving the expression of pro-cancer genes, which may regulate cell proliferation and migration. In most insect species, miR-100 clusters with miR-125 in the same primary transcript. These two miRNAs are involved in developmental timing in C. elegans and D. melanogaster (Rubio and Belles 2013). In the cockroach B. germanica, depletion of miR-100 with specific anti-miRNAs in the last instar nymph may reduce adult wing size (Rubio and Belles 2013). In the fruit fly D. melanogaster, miR-125 extends the lifespan by repressing chinmo in adult brains. Since their concentrations significantly change during rapa whelk development, the functions of these miRNAs in this species are of great interest. In R. venosa, miR-9-y was found at a low level (15.01 TPM) in competent larvae but rose to a high level (132.39 TPM) in postlarvae. It may prevent programmed cell death during wing development in Drosophila metamorphosis by repressing Drosophila LIM-only (Bejarano et al. 2010). This function may explain our finding that the miR-9 expression level was low in competent larvae. This developmental stage involves considerable amounts of programmed cell death because the velum must degenerate. Therefore, the levels of miR-9 must be high after metamorphosis in order to suppress further apoptosis.
As stated in previous studies (Song et al. 2016b,c), the mRNAs/ proteins involved in ingestion and digestion, cytoskeleton and cell adhesion, and apoptosis may also have important roles in driving metamorphosis. Twenty key miRNA-target pairs implicated in these processes were identified, and nine of them were further studied by real-time PCR. Both high-throughput sequencing and real-time PCR showed that the let-7-y expression level continuously increased during metamorphosis, whereas those of its target miRNAs like SARP-19, conotoxin, and exoglucanase continuously declined. The SARP-19 gene was expressed highly in the gastropod larval digestive gland and was sensitive to metamorphic cues (He et al. 2014). Conotoxin is a group of neurotoxic peptides isolated from Conus venom. High conotoxin levels in R. venosa pelagic larvae indicated that this life stage in R. venosa is homologous with that of Conus. In the rapa whelk, however, it degenerated after metamorphosis (Song et al. 2016d). Exoglucanase, an important digestive enzyme in whelk pelagic larvae, sharply decreased when the whelk become carnivorous after metamorphosis. The digestion-related genes were negatively regulated by let-7-y, which implies that let-7 participates in digestive system changes during whelk metamorphosis. The let-7 miRNA also participates in metamorphosis n in many other animals. In Drosophila, for example, the loss of let-7 and miR-125 may delay the terminal cell-cycle exit in the wing and the maturation of neuromuscular junctions (NMJs) in the adult abdominal muscles. The maturation rate of abdominal NMJs was governed by let-7 during metamorphosis by regulating the expression of the ab gene (Caygill and Johnston 2008). In the silkworm Bombyx mori, let-7 regulates molting and metamorphosis. Decreased let-7 expression in the silkworm could increase the expression of its target genes FTZ-F1 and Eip74EF (key regulatory factors in the ecdysone pathway) and cause developmental arrest during the larval-larval and larval-pupal transitions (Ling et al. 2014). The development-related miR-276 may inhibit apoptosis in shrimp hemocytes (Yang et al. 2012). In this study, miR-276 was found to be significantly downregulated in competent larvae, whereas the caspase-3 gene was upregulated. Therefore, miR-276 may regulate apoptosis by targeting the caspase-3 gene. Coexpression studies of key miRNA targets revealed their potential roles in whelk metamorphosis, but the mechanisms by which these miRNAs regulate this developmental process have not been fully explored yet.
In conclusion, the present study provides the first global view of the changes in miRNA that occur during rapa whelk metamorphosis. A total of 195 miRNAs were significantly differentially expressed and their target mRNAs were identified. These molecules are responsible for morphological and functional changes in organs. Some miRNAs involved in ingestion and digestion, cytoskeleton and cell adhesion, and apoptosis during metamorphosis are of great interest and were listed and validated by real-time PCR. These results will provide a basis for understanding the molecular mechanisms involved in the regulation of gastropod metamorphosis.

ACKNOWLEDGMENTS
This research was supported by the National Natural Science Foundation of China (NSFC, Grant No. 31572636), the earmarked Figure 4 Real-time PCR analysis for nine target miRNA-mRNA pairs. The appearance of the same letter over matching bars means that there is no significant difference between the stages (P . 0.05). Different letters over matching bars mean significant differences (P , 0.05). Values indicated are means 6 SE (N = 5).