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

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-value <= 0.01). With the criteria of fold change ≥1.2 and< 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 (). q p 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 (). 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). Figure 1

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< 0.05) between ESI and ESA crabs, and most of them (24 miRNAs, 77.42%) were downregulated in ESA individuals (and). p Figure 2 Table S1

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 (and). The miRNA-protein network was composed of 26 miRNAs and 30 target proteins. Among the targets, several proteins related to methyl farnesoate (MF) synthesis [Farnesoic acid O-methyltransferase (FAMeT)], oocyte meiotic maturation (Cyclin B), calcium signal transduction [Calmodulin (CaM)], and energy metabolism [AMP-activated protein kinase (AMPK), Phosphoglycerate mutase 2 (PGAM2), Triosephosphate isomerase (TPI), Phosphopyruvate hydratase (PPH), and Methylmalonyl-CoA epimerase (MCE)] were identified. Figure 3 Table 2

Using the dual luciferase reporter assay, we validated direct interaction of five miRNA-target pairs (miR-263a/FAMeT, miR-4171/FAMeT, miR-2b/CaM, miR-317/CaM, and miR-466f-3p/CaM) (). The results show that the relative luciferase activity (firefly luciferase activity/Renilla luciferase activity) was significantly reduced after pmirGLO-FAMeT-3′UTR was cotransfected with miR-263a or miR-4171 and pmirGLO-CaM-3′UTR was cotransfected with miR-2b, miR-317, or miR-466f-3p, which indicates the interaction between the miRNAs and their targets. Figure 4

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 () (miR-263a/:< 0.01; miR-4171/:< 0.01; miR-2b/:< 0.01; miR-317/:< 0.01), and there is no significant correlation between expression levels of miR-466f-3p and(> 0.05).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.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. Figure 5 FAMeT p FAMeT p CaM p CaM p CaM p FAMeT CaM

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 (). Consistent with the results of Vg transcript levels in previous studies (,), 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, 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 (). 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. ^{11 38 39 40} P. trituberculatus

Methyl farnesoate (MF), a sesquiterpene compound, plays a crucial role in reproduction regulation of decapods (). 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 (,,). The biosynthetic pathway of MF is similar to the general mevalonate pathway for acyclic isoprenoids, and the final step is catalyzed by FAMeT (–). 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 (). MF level is controlled by both anabolism and catabolism. Interestingly, a key enzyme in MF degradation, namely juvenile hormone esterase-like 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 (). 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. ^{41 1 42 43 44 46 46 47}

Calcium signals are demonstrated to be essential for oocyte maturation in mammals, fish, mollusks, and crustaceans (–). 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 maturation-promoting 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 (–). It is reported that ESA results in an increased level of CaM mRNA in the ovary of the red swamp crayfish, and knockdown of CaM expression by RNA interference suppresses the ovarian maturation induced by ablation (). 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 prawnfinds 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 (). 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. ^{48 50 51 53 54 11} Procambarus clarkia Penaeus monodon

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 theovary. 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 (,). 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 (,). 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. P. trituberculatus ^{55 56 57 58}

Accumulating evidence suggests that miRNAs are critical in ovarian development of crustacean (,). 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 onmiRNA target prediction and inverse correlation of protein and miRNA expression, we generated a network encompassing 26 miRNAs, 30 target proteins, and 71 potential functional associations (). 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, andand between miR-2b, miR-317, and. 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. ^{59 60} Figure 3 in silico FAMeT CaM

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. Femaleat 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. P. trituberculatus

During the period of the experiment, all the crabs were fed daily at 17:00 with live Manila clam, 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. Water salinity, pH, and ammonia were determined with a salinity refractometer (AIAGO, Japan), pH meter (WTW, Germany), and Hypobromite methods, respectively. Ruditapes philippinarum −1

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.

Raw LC-MS/MS data was converted into MGF format with the exported ProteoWizard tool (). Proteins were identified with Mascot version 2.3.02 (Matrix Science, UK) againstovarian transcriptome. The proteins containing at least one unique set of spectra were used for the following quantification analysis with IQuant software (). The proteomic data set was deposited in the iProX database under the accession number IPX0002268000. Avalue <= 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. ^{61 62} P. trituberculatus q

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.

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. 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 () and mirdeep2 (). 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. P. trituberculatus ^{63 64}

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 (). Differential expression analysis of the two libraries was performed using the DEGseq R package ().values were adjusted usingvalue (). The criteria of< 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 thegonad transcriptome data were used to predict putative targets of the identified miRNAs with miRanda software (). ^{65 65 66 67} P q q P. trituberculatus

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) usingI andI restriction sites. The plasmids were transformed into TOP10cells 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 (). Sac Sal E. coli 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 (< 0.05). p

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 Cmethod (). All the primers are shown in Supplementary Materials (). The correlation in expression between the miRNAs and their targets was determined using Spearman's correlation analysis. T ⁶⁸ Table S3