Integrative Proteomic and MicroRNA Analysis: Insights Into Mechanisms of Eyestalk Ablation-Induced Ovarian Maturation in the Swimming Crab Portunus trituberculatus

Eyestalk ablation is the most common method to induce ovarian maturation in decapod crustacean aquaculture, but it jeopardizes broodstock survival and larvae production. It is important to understand the molecular basis underlying the maturation triggered by ablation and thereby develop an alternative measure for maturation manipulation. In this study, we investigate alterations of ovarian proteome and miRNA profile after ablation in a commercially important marine crab Portunus trituberculatus. Quantitative proteomic analysis using iTRAQ reveals that 163 proteins are differentially expressed following ablation, and modulation of methyl farnesoate metabolism and activation of calcium signaling may play important roles in the ovarian maturation induced by ablation. miRNA expression profiling identifies 31 miRNAs that show statistically significant changes. Integration of miRNA and proteome expression data with miRNA target prediction algorithms generates a potential regulatory network consisting of 26 miRNAs and 30 proteins linked by 71 possible functional associations. The miRNA-protein network analysis suggests that miRNAs are involved in promoting ovarian maturation by controlling expression of proteins related to methyl farnesoate synthesis, calcium signals, and energy metabolism. Experimental validation and temporal expression analysis indicate multiple miRNAs can act synergistically to regulate expression of Farnesoic acid O-methyltransferase and Calmodulin. Our findings provide new insights for elucidating the mechanisms underlying eyestalk ablation-induced ovarian maturation and could be useful for devising an alternative technique for manipulating reproduction in P. trituberculatus and other decapods.


INTRODUCTION
The eyestalk, where the X-organ/sinus gland (XO/SG) complex is located, is an important neuroendocrine system in crustaceans, which is the major site of production and storage for a variety of important neuropeptides, such as the gonad-inhibiting hormone (GIH), the molt-inhibiting hormone (MIH), the pigmentdispersing hormone (PDH), and the red pigment-concentrating hormone (RPCH) (1,2). Among these neuropeptides, GIH is considered to play a key role in regulating female crustacean reproduction by repressing ovarian maturation and spawning (3). In decapod aquaculture practice, removal of GIH by eyestalk ablation (ESA) has been the most common method to induce ovarian maturation (4). However, in addition to GIH, secretion of other neuropeptides from the eyestalk is also impaired by this method, which results in significant hormonal imbalance and undesired impacts on not only the survival of broodstock, but also the quantity and quality of larvae produced (5). Hence, an alternative technique for maturation manipulation without ESA has been a long-term goal for the decapod culture industry.
To develop new techniques for triggering maturation, it is important to understand the molecular mechanisms of ablationinduced ovarian maturation. During the past decade, many efforts have been made to elucidate the molecular effects of ESA on ovarian development in decapods (6)(7)(8)(9)(10). With highthroughput transcriptome analysis, recent research provides a global picture of gene expression changes at the transcriptional level after ablation. Uawisetwathana et al. (11) find that eyestalk ablation activates several signal pathways important for ovarian development, such as the gonadotropin-releasing hormone (GnRH) signaling pathway and the progesterone-mediated oocyte maturation pathway, which promotes ovarian maturation in Penaeus monodon. Lee et al. (12) report that ESA can modulate gene expression in the hepatopancreas of Litopenaeus vannamei, upregulating genes related to metabolic processes as well as those associated with immunity and stress responses. These studies expand our knowledge on the transcriptional regulation that mediates ESA-induced maturation. Nevertheless, transcription change does not depict the changes of protein levels with complete accuracy because gene expression can also be modulated at the post-transcriptional and translational levels. To date, no information is available regarding gene expression regulation at post-transcriptional and translational levels following ESA.
MicroRNAs (miRNAs), a group of small non-coding RNAs, function as key post-transcriptional regulators of gene expression by complementary binding to 3 ′ untranslated regions of their target mRNA, leading to mRNA destabilization, or protein translation blockage (13,14). They are involved in many fundamental biological processes, such as development, metabolism, cell proliferation, and signal transduction (15)(16)(17). Recent evidence has revealed that miRNAs are critical in ovarian development and maturation in decapod crustaceans. In Eriocheir sinensis, miR-2, and miR-133 are demonstrated to regulate meiotic oocyte maturation by modulating cyclin B expression (18). A number of miRNAs have been found to exhibit ovary-biased expression in L. vannamei, Macrobrachium nipponense, Portunus trituberculatus, and Scylla paramamosain, and many of the predicted target genes for these miRNAs are crucial in regulating ovarian development, indicating the important functions of miRNA in ovarian development (19)(20)(21)(22). However, the roles of miRNA in ESA-induced ovarian maturation is still unknown.
Proteomic analysis is a tool that facilitates research on comprehensive protein expression profiling and identification of individual proteins involved in specific biological responses. It has been widely utilized in studying reproductive development in mammals, fish, and insects (23)(24)(25). In crustaceans, only a few proteomic studies have been conducted on ovarian development, and these studies all rely upon 2-D gel electrophoresis data (26)(27)(28), which cannot identify low abundant proteins, proteins with low or high molecular weights, and proteins that are excessively acidic or basic as well as hydrophobic proteins (29,30). The isobaric tags for relative and absolute quantitation (iTRAQ) is a mass spectrometry-based proteomics technique established in recent years and has overcome some of the drawbacks of 2-D gel electrophoresis; it is sensitive and can detect large proteins and low-abundant proteins with high throughput and low experimental error (31). Hence, iTRAQ should be an ideal approach to study protein expression changes after ESA.
The swimming crab P. trituberculatus (Crustacea: Decapoda: Brachyura) is an important fishery and aquaculture species, widely distributed in the estuary and coastal waters of Korea, Japan, China, and Southeast Asia (32). This species supports a large aquaculture industry in China, and its farming developed quickly during the last decade. In 2017, the production of P. trituberculatus reached 119,777 t (China Fishery Statistical Yearbook, 2018). However, due to rapid expansion of crab farming and lack of an effective technique for reproductive manipulation, high-quality seed production cannot fulfill the demand, which restricts further development of the industry (33,34). Therefore, it is crucial to understand the regulatory mechanisms of ovarian maturation and improve the means of artificial seed production. In this study, we investigate the alteration in proteome and miRNA transcriptome following ESA. The results can improve our understanding of molecular mechanisms underlying ESA-induced ovarian maturation and provide useful information for developing new reproduction manipulation techniques.

Protein Expression Altered by ESA
To identify proteins involved in ESA-induced ovarian maturation, iTRAQ was employed to assess ovarian protein expression changes at the fourth day after ESA. In total, 483,601 spectra were generated from the ovary of the eyestalk-intact (ESI) and ESA crabs. Based on the spectral data, 15,377 peptides and 3400 proteins were identified (cutoff: Mascot Percolator q-value <= 0.01). With the criteria of fold change ≥1.2 and p < 0.05, 163 proteins were identified to be differentially expressed (DE) between ESI and ESA crabs; a great majority (132 proteins, 80.98%) of the DE proteins showed upregulation in ESA individuals ( Table 1).    The functional category of the DE proteins was analyzed against the Gene Ontology (GO) database using three sets of ontologies: biological process, molecular function, and cellular component (Figure 1). The most abundant proteins in the biological process category were related to metabolic process (GO:0008152, 37 proteins), followed by cellular process (GO:0009987, 35 proteins). In the cellular component category, the proteins associated with cells (GO:0005623, 29 proteins) and cell parts (GO:0044464, 29 proteins) were dominant. In the molecular function category, binding (GO:0005488, 37 proteins) was the most prominent, followed by catalytic activity (GO:0003824, 34 proteins).

miRNA Expression Modulated by Eyestalk Ablation
To decipher the miRNAs implicated in ovarian maturation induced by ESA, the global ovarian miRNA expression after ESA was investigated using miRNA high-throughput sequencing. Relatively strict criteria were used to identify known and novel miRNAs. Only the small RNAs that aligned to the miRNAs in miRBase with no mismatch were classified as known miRNAs, and only those identified by both miREvo and mirdeep2 software were considered as novel miRNAs. In total, 184 unique miRNAs were identified after eliminating low-abundance miRNAs (reads ≤ 20), among which 100 were novel. Despite stringent criteria used for identifying miRNAs, there is possibility that some of these miRNAs may be degraded products of other RNAs due to lack of reference genome for miRNA prediction.
miRNA profile analysis showed that 31 miRNAs exhibited differential expression (under the criteria of fold changes ≥ 2.0 and p < 0.05) between ESI and ESA crabs, and most of them (24 miRNAs, 77.42%) were downregulated in ESA individuals (Figure 2 and Table S1).

In silico Correlation of Protein and miRNA and Experimental Validation
To uncover the potential crosstalk between specific miRNAs and proteins in response to ESA, we constructed a regulatory network based on miRNA target prediction algorithms and inverse correlation of protein and miRNAs expression (Figure 3 and Table 2  Phosphoglycerate mutase 2 (PGAM2), Triosephosphate isomerase (TPI), Phosphopyruvate hydratase (PPH), and Methylmalonyl-CoA epimerase (MCE)] were identified.
To further validate the correlation between the miRNAs and their targets and investigate their temporal express pattern after ESA, we analyzed the levels of the miRNAs and their target genes at different times following ESA. The results show that expression of miR-263a, miR-4171, miR-2b, and miR-317 negatively correlates with that of their target genes ( Figure 5) (miR-263a/FAMeT: p < 0.01; miR-4171/FAMeT: p < 0.01; miR-2b/CaM: p < 0.01; miR-317/CaM: p < 0.01), and there is no significant correlation between expression levels of miR-466f-3p and CaM (p > 0.05). FAMeT increased significantly from 72 h and maintained high expression until 168 h after ESA, whereas miR-263a and miR-4171 exhibited lower expression at those time points. CaM showed a significant upregulation from 48 to 168 h, and miR-2b and miR-317 were downregulated during these time periods. These results further confirm that miR-263a, miR-4171, miR-2b, and miR-317 are involved in regulating the expression of FAMeT and CaM after ESA.

DISCUSSION
Despite the undesirable effects, eyestalk ablation is currently the most effective method to induce ovarian maturation in the commercial hatchery of decapods. To date, the molecular mechanisms underlying the induced maturation have not been fully understood, which hinders development of new techniques for maturation manipulation. Previously, many studies have been focused on transcriptional regulation in response to ESA (35)(36)(37). Considering the fact that proteins and miRNAs form the cellular end point (proteins) of phenotypic features and their regulatory elements (miRNA), the present study concentrated on changes that occur in protein and miRNA expression. The results here present for the first time extensive ovarian proteome and miRNA profiling reflecting the alterations induced by ESA. Proteomic analysis using iTRAQ shows that 163 proteins are differentially expressed in the ovary of intact and ablated crabs, and among them, proteins regulating vitellogenesis, oocyte meiosis, and energy metabolism are identified, indicating their involvement in promoting ovarian maturation after ESA. miRNA transcriptome results show that 31 miRNAs exhibit differential expression after ESA. By integrating bioinformatic algorithms with expression data of miRNA and protein, we identify 30 targets of the DE miRNAs and define the potential roles of miRNA in maturation induced by the ablation.

Vitellogenesis
Vitellogenin (Vg), the precursor of the major egg yolk proteins, is regarded as a reliable marker of ovarian maturation in decapods as its expression is positively correlated to ovarian maturation levels (11). Consistent with the results of Vg transcript levels in previous studies (38,39), Vg protein level in the ovary was enhanced significantly after ESA in this study, confirming that ESA leads to increased ovarian vitellogenesis in the swimming crab. In crustaceans, various Vg synthesis sites have been reported in different species. For P. trituberculatus, vitellogenesis is demonstrated to occur both within the oocyte (endogenous vitellogenesis) and external to the oocyte (exogenous vitellogenesis), and exogenous Vg is the major source of Vg, which is synthesized in hepatopancreas, transported via hemolymph, and absorbed by oocytes or follicle cells through receptor-mediated endocytosis (40). In our study, an elevated level of Vg receptor (VgR) was observed in the ovary of the ablated crabs compared with that of the intact crabs. This result suggested that ESA can promote accumulation of exogenous vitellogenesis by upregulating VgR expression and facilitating uptake of Vg from hemolymph.

Methyl Farnesoate Metabolism
Methyl farnesoate (MF), a sesquiterpene compound, plays a crucial role in reproduction regulation of decapods (41). Previous studies demonstrate that there is a positive correlation between MF titer and ovarian maturation in decapods, and injection of MF stimulates oocyte growth and vitellogenesis (1,42,43). The biosynthetic pathway of MF is similar to the general mevalonate pathway for acyclic isoprenoids, and the final step is catalyzed by FAMeT (44)(45)(46). In this study, a significant upregulation of FAMeT in ovary was observed after ESA, indicating that ESA resulted in an increased rate of MF synthesis. This result is in accordance with previous studies that find the MF levels are negatively regulated by eyestalk neuropeptides, and removal of the eyestalk leads to elevated levels of MF (46). MF level is controlled by both anabolism and catabolism. Interestingly, a key enzyme in MF degradation, namely juvenile hormone esteraselike carboxylesterase (JHE-like CXE), was also found upregulated in ESA crabs, indicating an increase in MF catabolism in addition to MF anabolism. It is reported that a high level of MF has a deleterious effect on oocyte growth and development, resulting in oocyte degeneration (47). The increased MF catabolism in the ovary after ESA may represent a protective mechanism to regulate MF levels, avoiding the detrimental effects of excessive MF.

Calcium Signals
Calcium signals are demonstrated to be essential for oocyte maturation in mammals, fish, mollusks, and crustaceans (48)(49)(50). Elevated concentration of intracellular calcium in the oocyte can activate CaM, a molecular switch to regulate the network of calcium signaling and, in turn, trigger a cascade of downstream signaling, leading to activation of maturationpromoting factor (MPF), which is a heterodimeric complex composed of Cell division cycle 2 (Cdc 2) and Cyclin B and responsible for triggering G2/M-phase transition in oocytes (51)(52)(53). It is reported that ESA results in an increased level of CaM mRNA in the ovary of the red swamp crayfish Procambarus clarkia, and knockdown of CaM expression by RNA interference suppresses the ovarian maturation induced by ablation (54). In this study, CaM, as well as Cdc 2 and Cyclin B, were upregulated after ESA. These results together with previous findings suggest that ESA can stimulate meiotic maturation through calcium signaling. In mammals, GnRH signaling can trigger calcium signaling, thereby inducing oocyte maturation. Similarly, a recent study in the black tiger prawn Penaeus monodon finds that genes in the GnRH pathway exhibit an early induction pattern after ESA and speculates that GnRH signaling leads to the activation of calcium signaling (11). However, upregulation of the proteins in the GnRH signaling pathway was not observed in this study. It is possible that activation of GnRH signaling is transient at the onset of oocyte maturation after ESA. Further study is required to confirm the upstream signals triggering calcium-signaling pathways after ESA. Metabolism KEGG analysis shows that the proteins in metabolic pathways account for the largest proportion of DE proteins (22.70%), and they are mapped to a wide range of pathways related to metabolism, such as carbohydrate, amino acid, energy, lipid, and nucleotide metabolism, indicating that ESA has an extensive impact on the metabolism in the P. trituberculatus ovary.
It is noteworthy that a significant upregulation of AMPK was observed after ESA. As a key regulator that maintains cellular energy homeostasis, AMPK reprograms cellular metabolism from anabolism to catabolism by controlling lipid and glucose metabolism when activated by falling energy status (55,56). In accordance with AMPK upregulation, several important enzymes in glycolysis (PGAM2, TPI, and PPH), tricarboxylic  acid cycle (ACO), and fatty acid catabolism (MCE) were also found upregulated following ablation. Previous studies in other decapods show that ESA can promote the transportation of glucose and lipid from the hepatopancreas to the ovary (57,58). Those findings and our results together indicate ESA results in an increased energy requirement for the accelerated ovarian maturation. Hence, in aquaculture practice, the feeds for the broodstocks should be adjusted to meet the high energy demand after ablation.

The Roles of miRNA
Accumulating evidence suggests that miRNAs are critical in ovarian development of crustacean (59,60). However, their functions in ESA-induced ovarian maturation are still unclear.
In this study, 31 miRNAs were differentially expressed in ESA and ESI crabs. To uncover the functions of these miRNA during accelerated ovarian maturation induced by ESA, we integrated miRNA transcriptome with proteome data sets. Based on in silico miRNA target prediction and inverse correlation Above the histograms, nucleotide sequences of miRNAs and targets wild-type and mutant 3 ′ UTR regions are shown. The histograms present relative luciferase activity in HEK293T cells cotransfected with miRNA mimics or miRNA negative control, and pmirGLO vector containing wild-type or mutant 3 ′ UTR regions of miRNA targets. The relative luciferase activity was calculated by normalizing firefly luciferase activity to Renilla luciferase activity. Different letters on the bars indicate significant differences based on one-way ANOVA followed by Tukey's HSD multiple comparison test (p < 0.05). of protein and miRNA expression, we generated a network encompassing 26 miRNAs, 30 target proteins, and 71 potential functional associations (Figure 3). Among the targets, a number of proteins mentioned above, including FAMeT, Cyclin B, CaM, AMPK, TPI, and PPH, were identified, suggesting that miRNAs may play important roles in ESA-mediating ovarian maturation through regulating proteins associated with MF synthesis, calcium signals, and energy metabolism. Within the miRNA-protein network, most miRNAs have multiple targets, and inversely, many proteins are targeted by several miRNAs, implying a complex post-transcriptional regulation between miRNAs and proteins in response to ESA. To experimentally validate the synergistic target regulation of the miRNAs, we tested FAMeT regulation by miR-263a and miR-4171, and CaM regulation by miR-2b, miR-317, and miR-466f-3p, using the dual luciferase reporter assay. The results indicate that both miR-263a and miR-4171 directly regulate FAMeT expression, and miR-2b, miR-317, and miR-466f-3p repress expression of CaM. In addition, we investigated the expression pattern of the miRNAs and their target genes after ESA, and the results show significant negative correlation between miR-263a, miR-4171, and FAMeT and between miR-2b, miR-317, and CaM. Taken together, these data demonstrate that multiple miRNAs could function cooperatively to regulate FAMeT and CaM expression, and thereby participate in inducing ovarian maturation after ESA.

CONCLUSIONS
In summary, we characterize the changes in ovarian proteome and miRNA transcriptome in P. trituberculatus after ESA.
The results suggest that alterations in MF metabolism and calcium signaling are crucial for promoting ovarian maturation, and ESA results in a higher energy production to meet the increased energy demand in the ovary. Furthermore, our study reveals a miRNA-mediated mechanism for inducing ovarian maturation. These findings improve our understanding of the molecular mechanisms of ESA-induced ovarian maturation and lay the foundation for developing an alternative technique for maturation manipulation without ESA in decapod crustaceans.

Animal Collection and Eyestalk Ablation
All the experimental procedures involving the handling and treatment of the crabs used in this study were approved by the Institutional Animal Care and Use Committee of Yellow Sea Fisheries Research Institute prior to initiation of experiments. Female P. trituberculatus at 6-month age (198.65 ± 23.77 g) were collected from Haifeng Company, Weifang, China. The crabs were acclimated at laboratory conditions (temperature, 21 • -23 • C; salinity, 30-31) for 3 weeks. Then, three crabs were bilaterally eyestalk-ablated, and wounds were cauterized to minimize loss of hemolymph. Based on the results of our preliminary experiment, 96 h post ESA was selected as sample time; three eyestalk-ablated individuals (ESA) and three individuals with intact eyestalks (ESI) were placed in an ice bath until anesthetized, and the ovaries were collected and immediately frozen in liquid nitrogen and stored at −80 • C. During the period of the experiment, all the crabs were fed daily at 17:00 with live Manila clam Ruditapes philippinarum, and the feces and leftover feed were removed prior to feeding. Aeration was provided continuously, and the photoperiod was 12 h light:12 h dark. Seawater was filtered using a sand filter, and one third to one half of the rearing water was exchanged using fresh equi-temperature seawater. Water pH was around 7.5, and ammonia was <0.23 mg L −1 . Water salinity, pH, and ammonia were determined with a salinity refractometer (AIAGO, Japan), pH meter (WTW, Germany), and Hypobromite methods, respectively.
Protein Extraction, Digestion, and iTRAQ Labeling iTRAQ analysis was performed at Beijing Genomics Institute (BGI, Shenzhen, China). The ovary of each crab from the ESI and ESA groups were disrupted in lysis buffer with enzyme inhibitors by TissueLyser (Qiagen, USA). The mixtures were centrifuged at 25,000 g for 20 min, and the supernatant was carefully removed and mixed with 5 volume of cold acetone and stored at −20 • C for 2 h prior to centrifuging. The pellets were dissolved with lysis buffer, and 10 mM dithiothreitol (DTT) was added and maintained at 56 • C for 1 h to reduce the disulfide bond of peptides. Then, 55 mM IAM was added to the solution and kept in a dark room for 45 min. After adding 5 volume of chilled acetone into the solution and kept at −20 • C for 2 h, the solution was centrifuged again, and the pellet was dissolved with lysis buffer to get a protein solution. The protein concentration was determined using the Bradford method.
The protein solutions (100 µg) from each sample were digested with Trypsin gold (Promega, USA). After digestion, the peptides were vacuum centrifuged to dryness and dissolved with 0.5 M TEAB. The iTRAQ labeling of peptides was performed on the ESI and ESA groups using the iTRAQ Reagent 8-plex kit (Applied Biosystems, USA) according to the manufacturer's protocol.

Protein Identification, Quantification, and Functional Analysis
Raw LC-MS/MS data was converted into MGF format with the exported ProteoWizard tool (61). Proteins were identified with Mascot version 2.3.02 (Matrix Science, UK) against P. trituberculatus ovarian transcriptome. The proteins containing at least one unique set of spectra were used for the following quantification analysis with IQuant software (62). The proteomic data set was deposited in the iProX database under the accession number IPX0002268000. A q value <= 0.05 and foldchange >= 1.2 were set as the threshold for differentially expressed proteins (DEPs). The GO and KEGG databases were used to classify and group the DEPs.

RNA Extraction, Library Construction, and miRNA Sequencing
Total RNA was extracted from the ovary of the same crabs used for proteomic analysis with TRIzol reagent (Invitrogen, USA) and purified with a mirVana miRNA Isolation Kit following the manufacturer's protocol (Ambion, USA). Total RNA degradation and contamination was accessed on 1% agarose gels. Quantity and integrity of the RNA samples were determined using a Nano Photometer spectrophotometer (Implen, USA) and Bioanalyzer 2100 system (Agilent Technologies, USA). The RNA concentration was measured using the Qubit RNA Assay Kit in the Qubit 2.0 Flurometer (Life Technologies, USA).
For the ESI or ESA small RNA library, an equal amount of RNA samples from each replicate were pooled together. The sequencing libraries were generated using NEBNext Multiplex Small RNA Library Prep Set for Illumina (New England Biolabs, USA) according to the manufacturer's protocol, and index codes were added to attribute sequences to each sample. Library quality was assessed on the Agilent Bioanalyzer 2100 system (Agilent Technologies, USA) using DNA High Sensitivity Chips. The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq SR Cluster Kit v3-cBot-HS (Illumina, USA) according to the manufacturer's instructions. After cluster generation, the small RNA libraries were sequenced on an Illumina Hiseq 2500 platform at Novogene Company, Beijing, China. The small RNA sequencing data set was deposited in NCBI Sequence Read Archive (SRA) under the accession numbers PRJNA639350.

miRNA Data Analysis
After Illumina sequencing, clean reads were obtained by removing low-quality sequences, adapter-contaminated tags, and reads with poly N (where N represents unknown base). All the clean reads were searched against GenBank and Rfam databases to exclude known non-coding RNAs, including rRNAs, tRNAs, snRNAs, and snoRNAs. Any reads encoding proteins were also removed by blasting against the reference unigenes derived from the gonadal transcriptome data set of P. trituberculatus. Then, the remaining sequences were searched against the miRNAs from all the animals in miRBase to identify known miRNAs. Only the sequences that aligned to the miRNAs in miRBase without any mismatch were considered as known miRNAs. To predict novel miRNAs, the remainder unannotated small RNA sequences were analyzed with an integrated combination of miREvo (63) and mirdeep2 (64). The software identified novel miRNAs by exploring the secondary structure, the Dicer cleavage site, and the minimum free energy of the small RNA tags unannotated in the former steps, and the threshold score set as ≥ 50. The miRNAs predicted by the software were considered as novel miRNAs.
To analyze the expression profiles of the miRNAs in the ovary of ESA and ESI crabs, the read counts of miRNAs were normalized into TPM (transcript per million) through the normalization formula: Normalized expression = (Mapped reads/Total reads) ×106 (65). Differential expression analysis of the two libraries was performed using the DEGseq R package (65). P values were adjusted using q value (66). The criteria of q < 0.01 and |log2(fold-change)| > 1 was set as the threshold for defining statistically different expression. Because genome information for the swimming crab is not available, 3 ′ UTR sequences extracted from the P. trituberculatus gonad transcriptome data were used to predict putative targets of the identified miRNAs with miRanda software (67).

miRNA Target Validation
Five selected miRNA-target pairs (miR-263a/FAMeT, miR-4171/FAMeT, miR-2b/CaM, miR-317/CaM, and miR-466f-3p/CaM) were validated with dual luciferase reporter assay. The 3 ′ UTR sequences of FAMeT and CaM containing wild-type and mutant miRNA binding sites were artificially synthesized by General Biol Co. Ltd. (Hefei, China) and cloned into pmirGLO dual-luciferase reporter vector (Promega, USA) using SacI and SalI restriction sites. The plasmids were transformed into TOP10 E. coli cells and purified with a TIANprep Mini Plasmid Kit (Tiangen, China), and all the insertions were certified by DNA sequencing. The synthesized 3 ′ UTR sequences of FAMeT and CaM are shown in the supplementary materials (Table S2).
For the luciferase reporter assay, HEK293T cells were seeded in 24-well plates and transfected with 50 nM miRNA mimics or scrambled miRNA, and 1 µg luciferase reporter plasmid pmirGLO-wild type or pmirGLO-Mutant using Exfect 2000 Transfection Reagent (Vazyme, China). At 48 h after transfection, firefly and Renilla luciferase activities were determined using the dual luciferase reporter assay system (Promega, USA). The firefly luciferase signal was normalized to the Renilla luciferase signal. All experiments were performed in three replicates. The normalized firefly luciferase activity was compared between different groups using Tukey's HSD multiple comparison test (p < 0.05).
To further validate the expression correlation between the five selected miRNA target pairs, we analyzed levels of the miRNAs and their target at different time (0, 24, 48, 72, 96, and 168 h) after ESA. Because there is no proper antibody available, we measured mRNA levels of the targets instead. Total RNA of the samples was extracted using the RNAprep Pure Kit (Tiangen, China), and reverse transcription and RT-PCR of the target genes were performed with the FastKing RT Super Mix (Tiangen, China) and SuperReal PreMix Plus Kit (Tiangen, China), respectively, according to the manufacturer's protocol. The RT-PCR was programmed at 95 • C for 15 min, followed by 40 cycles of 95 • C for 10 s, 58 • C for 20 s, and 72 • C for 30 s. miRNA was isolated using the miRcute miRNA Isolation Kit (Tiangen, China) following the manufacturer's instruction, and reverse transcription and RT-PCR of the miRNAs were performed using miRcute miRNA First-strand cDNA Synthesis Kit (Tiangen, China) and miRcute Plus miRNA qPCR Kit (Tiangen, China). The RT-PCR was programmed at 95 • C for 15 min, followed by 40 cycles of 94 • C for 20 s and 60 • C for 34 s. Three biological replicates were measured at each time point, and each measurement was performed in triplicate. Relative levels of miRNAs and mRNAs were normalized to the U6 snRNA and β-actin, respectively, in each sample using the comparative C T method (68). All the primers are shown in Supplementary Materials ( Table S3). The correlation in expression between the miRNAs and their targets was determined using Spearman's correlation analysis.

DATA AVAILABILITY STATEMENT
The small RNA sequencing dataset was deposited in NCBI sequence read archive (SRA) under the accession numbers PRJNA639350.

ETHICS STATEMENT
The animal study was reviewed and approved by the Institutional Animal Care and Use Committee of Yellow Sea Fisheries Research Institute.

AUTHOR CONTRIBUTIONS
XM and PL: conceptualization, resources, supervision, project administration, and funding acquisition. MZ, XM and PL: methodology and investigation. XM, MZ, and JL: software, validation, formal analysis, data curation, and original draft. XM, MZ, and PL: writing, review, and editing.