All animal experiments in this study were conducted in accordance with the relevant national and international guidelines. Our project was approved by the East China Sea Fisheries Research Institute. In China, catching wild mud crabs from seawater does not require specific permits. Our study did not involve endangered or protected species.

Five healthy female S. paramamosain weighing 300 ± 30 g after reproductive molting were captured from the cultured population in the Ninghai Center of East China Sea Fisheries Research Institute. Once dissected, the hepatopancreas tissue samples of each crab were cut into 15 pieces, and were taken for in vitro experiments. Each female provided 3 similar pieces for each treatment group (random block design).

Each piece of hepatopancreas was placed in a well of a 12-well sterile culture plate (Corning Co. Ltd.^®). Each well was previously filled with 2 mL of culture medium DMEM (Gibco) at 1g/L glucose with enhanced amino acids and vitamins, without proteins and lipids. A hormone concentration of 1 μM, just able to stimulate JH signaling pathway genes in the hepatopancreas in preliminary trials, was added. The treatments were as follows: control (n=15); MF at 1 μM (n=15); FA at 1 μM (n=15); JH III at 1 μM (n=15); Met at 1 μM (n=15). MF and FA were purchased from Echelon Biosciences, Salt Lake City, UT, USA, and JH III and Met from Sigma-Aldrich. The plates were incubated for 6 h in a shaker, under room temperature. At the end of the incubation, all the hepatopancreas pieces were suddenly frozen in liquid nitrogen before storage at -80°C until metabolomics analysis. Mixtures 1-4 were prepared as quality control (QC) to test the repeatability of analysis. The triplicate samples from the same crab hepatopancreas in each treatment group were homogenized as one sample for further analysis.

Hydrophilic compounds: Samples were thawed on ice. 50 mg of each sample was homogenized with 1000 μL ice-cold methanol/water (70%, v/v). Cold steel balls were used to the mix and homogenize each sample at 30 Hz for 3 min. The mixture was whirled for 1 min, and then centrifuged at 12,000 rpm, 4°C, for 10 min. The supernatant was collected for LC-MS/MS analysis.

Hydrophobic compounds: Samples were thawed on ice. 50 mg of each sample was homogenized with 1 mL mixture (including methanol, MTBE and internal standard mixture) using steel ball. The mixture for another 2 min. 500 μL water was added and whirled for 1 min, and then centrifuged at 12,000 rpm, 4°C, for 10 min. 500 μL supernatant was extracted and concentrated. The powder was dissolved with 100 μL mobile phase B, and the dissolving solution was put into the sample bottle for LC-MS/MS analysis.

Analyst 1.6.3 was used to analyze the metabolomics data. Based on the self-built target database MWDB (metware database), qualitative analysis was carried out according to the retention time of substances, ion pair information and secondary spectrum data. Metabolite quantification was completed by multiple reaction monitoring (MRM) analysis of triple quadrupole mass spectrometry. In MRM mode, the first four-stage rod first screened the precursor ions of the target substance and eliminated the precursor ions corresponding to other substances to preliminarily eliminate interference. The precursor ions were broken after collision induced ionization in the second four-stage rod, and broke according to the structural characteristics of the material to form a series of unique fragment ions. The fragment ions were filtered through the third four-stage rod to select a typical characteristic fragment ion, without interference of the untargeted ions, and the quantification was more accurate with better repeatability. After obtaining the data of different samples, MultiQuant software was used to extract ion chromatographic peaks of all metabolites, and to integrate under the peak area, and the chromatographic peaks of the same metabolite in different samples were integrated and corrected.

The raw data was first unit variance scaled to prevent bias from large variance of the higher values. Simca 14.1 was used to perform Principal Component Analysis (PCA) for clustering trends, with no outliers excluded in our data. To improve group separation, Orthogonal Projections to Latent Structures Discriminant Analysis (OPLS-DA) was also applied, providing better separations between treatment groups (44). Model validity of the PCA and OPLS-DA models were checked via 7-fold cross-validation, which generated R2 (R2X and R2Y) and Q2 values to describe the goodness of fit and the predictive ability, respectively. Potential overfitting of the OPLS-DA model was checked by performing permutation tests. Differential metabolites were selected according to their variable importance in the projection (VIP) values from the OPLS-DA model, plus the fold change (FC) values calculated. The threshold was set at FC ≥2 and ≤0.5, and VIP≥1 to generate a list of significant discriminant metabolites between different treatments. The online database of Kyoto Encyclopedia of Genes and Genomes (KEGG; www.genome.jp/kegg/) was used for metabolite identification, biological role explanation, and pathway construction (45), and most of the affected genes and pathways are lipid metabolism-related.

Key regulators with important roles in the production of significant different metabolites in lipid metabolism pathways were screened and the expression patterns for the encoding genes were investigated in all groups. Total RNA was was extracted from hepatopancreas samples using Total RNA Extraction Reagent (EZ-10 DNAaway RNA Extraction Kit, BBI Life Scinece), after which the quantity and quality of total RNA were assessed by NanoDrop 2000 spectrophotometer. Then, 0.2ugRNA was converted into cDNA (complementary DNA) by ReverTra Ace qPCR RT Kit (TOYOBO). Quantitative real-time polymerase chain reaction (qPCR) was performed using PerfectStart Green qPCR SuperMix (TransGen Biotech). The amplification system was performed in a 10 μL reaction volume, containing 5 μL of 2×PerfectStart Green qPCR SuperMix, 0.2 μL of forward primer (10 μM), 0.2 μL of reverse primer (10 μM), 2 μL of cDNA and 2.6 μL nuclease-free water. The reaction program for qPCR was performed as follows: initial denaturation at 50°C for 2min and then 95°C for 10 min, 40 cycles at 95°C for 5s, 60°C for 30s, followed by a continuous increase to 95°C for melting curve analysis. The gene encoding 18S rRNA was used as reference. The selected genes were those encoding hepatocyte nuclear factor 4 (HNF4), carnitine palmitoyl transferase (CPT), very-long-chain acyl-CoA dehydrogenase (VLCAD), palmitoyl protein thioesterase (PPT) and 3-ketoacyl-CoA thiolase (KAT). Table 1 shows the primers used according to genome study (46). Three replicates were prepared for each sample. mRNA expression levels were determined using the 2^−ΔΔCT method (47). For each gene, the relative quantity of mRNA in arbitrary unit was interpreted as FC in comparison with the control group.

For the FCs in the raw intensities of significant metabolites and the mRNA relative expression levels of the selected genes, due to the non-parametric distribution, independent sample Kruskal Wallis tests were performed with IBM SPSS 26 to compare the differences among groups at P ≤ 0.05 as plotted in the box-and-whiskers.

In total, 1258 raw-data peaks were identified in hepatopancreas samples by four different hormonal treatments (FA, JH III, MF, Met) and the control group, including 474 substances from metabolomics and 784 substances from lipidomics. Representative UPLC-MS/MS total ion chromatograms (TICs) from the QC groups were presented in Figure A1 . After integral correction, all 1258 peaks were retained ( Figure A2 ). Afterwards, Z-score standarization was applied due to fairly different ranges of the original values.

The unsupervised PCA of the hepatopancreas metabolic profiles is presented in Figure A3 . The validity of the PCA model was evaluated in terms of R2X and Q2 values, 71.3% and 54.5% respectively, with values >0.5 indicative of high accuracy. However, the different treatments were not as clearly separated as different crab samples, indicating very similar impacts of the four hormonal subjects. Another explanation is that better separation should be obtained from a larger sample size to reduce irrelevant noise signals. To further differentiate between the treatment groups, as shown by the supervised OPLS‐DA scores plot ( Figure 2 ), the treatments were all clearly separated from each other on the X axis, suggesting metabolic perturbation caused by different hormonal applications. Two trends were identified from the X axis, with grouping of control and FA, as opposed to JH III, MF and Met, as separated by the Y axis. After 200 permutations for each treatment group, OPLS-DA was verified with all the intercepts of R2 on Y axis less than 0.3, and all the intercepts of Q2 on Y axis less than 0, indicating no overfitting. The three key parameters R2X, R2Y and Q2 were 0.80, 0.24 and -0.138 respectively. The R2X showed great accuracy of the model on the X axis, yet the Q2 value implied high similarity between hormonal treatments which coincides with the PCA result, indicative of their functional similarities. Both PCA and OPLS-DA results show that the QC samples are tightly gathered together showing good test quality and reliability of data, which confirms that the differences observed were more likely to result from genuine metabolite changes rather than technical errors.

One-to-one group comparison OPLS‐DA plots are shown in Figure A4 for differential metabolites. The number of significantly different metabolites (FC ≥2 and ≤0.5, VIP≥1) in different group comparisons is shown in Table 2 . The differences in metabolites were further analyzed between FA and control, JH III and control, MF and control, and Met and control, in order to clearly identify the effects of each hormonal treatment. The most differential metabolites were characterized between MF and control, of which 21 were upregulated and 2 were downregulated. This is followed by Met and control, with 18 compounds differentially expressed: 16 increased and 2 decreased. JH III and control comparison reported 16 differentially metabolized biomarkers with 13 increased significantly and 3 decreased significantly. Between FA and control, 12 metabolites differed significantly, of which 10 upregulated and 2 downregulated. The differential metabolites in the four comparison groups were ranked according to their FC value ( Tables A1 – A4 ) and visualized in volcano plots according to their VIP and FC values in Figures 3A – D . Venn diagram is also presented to filter common and unique metabolites with their up/down changing trends among the four comparison groups in Figure 3E . There were 2 metabolites common in all four treatments when compared with control: 3-N-methyl-L-histidine and D-arabitol; 8 metabolites common in JH III, MF and Met treatments compared with control: L-homocystine, lyso-phosphatidylcholine (LPC) 18:0 and 18:3, lyso-phosphatidyl ethanolamine (LPE) 18:0, lyso-phosphatidyl glycerol (LPG) 18:1, platelet activating factor PAF C-16 (1-O-hexadecyl-2-O-acetyl-sn-glyceryl-3-phosphorylcholine), sphinganine and tryptamine. Between FA and JH III, 1 common metabolite, (±)18-HEPE (EPA-metabolite 18-hydroxyeicosapentaenoic acid), differed significantly compared with control. Between MF—control and Met—control, LPC 15:0, 16:0 and 17:0, LPE 16:0, and triglyceride (TG) (18:0/20:0/22:0) are the five common differently expressed metabolites. Of the common significant metabolites, only L-homocystine was downregulated, while the others were all upregulated. The unique differently expressed metabolites for each treatment compared with control were as follows: allysine, cysteine glutathione disulfide, LPC 16:0, 20:3, 20:4, 22:4 and O-22:4, LPE P-18:0 and phosphatidyl ethanolamine (PE) P-18:1/18:0 only for FA application; 9-hydroperoxyoctadecadienoic acid (9-HPODE), free fatty acid (FFA) (20:4), O-phospho-L-serine, oleamide and oleate only for JH III application; 18-hydroxycorticosterone, 2-(dimethylamino)guanosine, lysophosphatidic acid (LPA) 16:0, LPE (16:1(9Z)), L-phenylalanine-L-phenylalanine, TG(16:0/20:0/22:0), TG(16:0/22:0/22:0) and TG(20:2/20:1/22:0) only for MF application; L-tartaric acid, N-lactoyl-phenylalanine and tryptophyl-glutamic acid only for Met application. Of the 41 significant metabolites in Venn diagram, 31 were lipids and involved in lipid metabolic pathways. Therefore, further analysis was focused on the predominant impacts of our hormones on lipid accumulation and related mechanisms.

The absolute quantifications and percentages of the overall lipids in the hepatopancreas of the mud crab were categorized in a pie chart in Figure 4A , with more than half of the lipid accumulated being TG(64%), followed by PC (9.65%) and LPC (5.27%). Altogether four types of TGs were identified as significant metabolites and their FC values were shown in Figure 4B . PC was not found in significant changing metabolites. LPC and other lysophospholipids accounted for more than half of the significant lipids (16 of the 31 total significant lipids), and their FC patterns were shown in Figure 4C . Thus, the major changing lipids we identified were TGs and lysophospholipids.

The relative mRNA expression of genes related to lipid metabolism in hepatopancreas of the mud crab were selected according to previous study (48) as presented in Figure 5 , in terms of both lipid degradation (A) (B) (C) (D) and lipid synthesis (E). The results represented the FCs of the expression of the target genes in all treatments against control. Three of the four lipid degradation genes (hnf4, cpt, ppt) were significantly upregulated by FA not JH III, MF or Met, opposite to the trend of TG accumulation, whereas the lipid synthesis gene (kat) was not stimulated by any treatment.

Several hormones have been implicated in the regulation/control of ovarian maturation in crustaceans including JH III, MF and FA (27). Our study focused on metabolomic and lipidomic changes among different juvenoids with potential MF activity, with overall 474 substances from metabolomics and 784 substances from lipidomics detected. An up-to-date study conducted using MS-based lipidomic approach in hepatopancreas for mud crab only detected a total of 390 lipid molecular species belonging to 22 lipid classes mainly glycerophospholipids, fatty acyls, glycerolipids and sphingolipids (49). Other study using untargeted metabolome analysis of hepatopancreas of the swimming crab Portunus trituberculatus detected a total of 628 metabolites by LC-MS (50). In hepatopancreas of Portunus trituberculatus, untargeted lipidomics identified a total of 1741 different lipid molecules from 32 lipid subclasses consisting of TGs, phosphatidylcholines, phosphatidylethanolamines and other lipid classes (51), which is consistent with the lipid categories we identified in Figure 4A . The enriched pathways in Table 2 were typical crab hepatopancreas pathways similar to a comparative proteomics analysis of the hepatopancreas at key ovarian stages in Eriocheir sinensis, with protein digestion, fatty acid metabolism, prostaglandin metabolism, lipid digestion/transportation and suchlike pathways differentially expressed (52).

Among terpenoid hormones, MF might be the principal effector in S. paramamosain ovarian maturation in our experiment with the most metabolites expressed and pathways activated compared with the control group ( Table 2 ). The OPLS‐DA clustering in Figure 2 indicated FA had less effect than JH III, MF and Met. This is congruent with the trend from the Venn diagram in Figure 3E , where altogether 8 metabolites were common in JH III, MF and Met treatments compared with control apart from the 2 common significant different metabolites for all four hormonal treatments, yet there is merely 1, 0 and 0 common metabolite between FA and JH III, MF, Met respectively, when compared with control. FA showed the largest number of unique significant different metabolites upregulated, indicative of a clearly differential role of this substance as opposed to JH III, MF and Met. Structurally, Figure 1 presented the unique carbonate acid group without esterification in FA, different from all other three sesquiterpenoids, denoting different specificities to different binding sites. Juvenile hormone acid methyltransferase (JHAMT) and farnesoic acid methyltransferase (FAMeT) are both involved in the MF biosynthesis from precursor FA in crustaceans. However, FAMeT activity in Metapenaeus ensis (53) and Cancer pagurus (54), JHAMT and FAMeT activity in Portunus trituberculatus (55), as requirements for the conversion of FA to MF, was solely detected in MO instead of hepatopancreas. But in Macrobrachium rosenbergii, FAMeT mRNA was highly expressed in the hepatopancreas (56). It is likely that FA is only produced and transformed to MF in MO not hepatopancreas in S. paramamosain (57). Therefore, a possible explanation is that the converting enzymes might be absent in the hepatopancreas, or inactive at the time of our hormonal addition, blocking downstream MF biosynthesis pathways from FA. This could be supported by the fact from KEGG annotation, that FA induced no other pathway apart from the only common pathway for all four hormones compared with control, pentose and glucuronate interconversions. Another explanation is the different effective concentrations. In the red crab Charybdis feriatus, distinct hormonal effects in vitro indicated the stimulatory effect of FA is more potent at low concentrations than MF and JH-III with regard to the vitellogenin gene upregulation (27). However, we only used one concentration for our experiments. Besides, MF shows inhibiting role in vitellogenesis in female crabs at middle and late vitellogenic stages, since the MF level is controlled by both anabolism and catabolism (58), whereas the role of FA reported is only for vitellogenic stimulation. The degradation pathway of MF may not be activated when there is only FA application, which also contributes to their differential effects.

The difference in metabolome and lipidome patterns following treatment with JH III, MF, FA and Met suggests that these hormones may be distinct in their effects during reproduction despite their close molecular similarity. The activity of JH III and Met should be more similar with MF activity, different from that of FA.

Figure 3E showed the significantly differed metabolites with the majority belonging to lipid metabolites, which mostly followed the grouping of control and FA versus JH III, MF and Met as denoted by shared letters in Figures 4B, C . Indeed, during the processes of reproduction, abundant lipids are accumulated and deposited in the hepatopancreas of crustaceans (59), either for energy supply or transferred to ovary to stimulate maturation, as reported in S. paramamosain (41). The diverse lipid categories in Figure 4A presented a large percentage at 64% in mud crab hepatopancreas was TG, followed by PC (9.65%) and LPC (5.27%). The significantly changed lipids we identified were also detected by Wang et al. in the hepatopancreas of the mud crab (49). Also, lipidomic analysis identified in Daphnia magna the changed levels of TGs, glycerophospholipids and cholesterol in females exposed to juvenoids (60). Among the significant lipid metabolites, TG and lysophospholipids were the major changing lipids. Figure 4B showed the four significantly altered TGs. The fatty acid compositions in TGs were the final reflection of fatty acid accumulation in crab hepatopancreas, with JH III, MF and Met significantly upregulated the types of TGs containing mostly long-chain saturated fatty acids such as 16:0, 18:0, 20:0 and 22:0 for energy storage in preparation for ovarian maturation. This is also supported by the fact that storage lipids were significantly enhanced by MF in Daphnia magna (61). However, FA increased TG accumulation to a much lesser extent.

LPA 16:0; LPC 15:0, 18:0, 18:3; LPE 16:1, 18:0; and LPG 18:1 also showed similar patterns under different treatments, with significant elevations in JH III, MF and Met rather than the FA and control group in Figure 4C , except for LPC 16:0 with a significant upregulation in FA group rather than JH III, MF and Met. Lysophospholipids are transferred from glycerophospholipids with the substitution of fatty acid chain in glycerophospholipids with hydroxyl group. Lysophospholipids are biologically active in a range of important cellular signaling pathways, functioning as ligands for G-protein-coupled receptors (62). In vitro and in vivo studies from the past decades have demonstrated or suggested the physiological functions of lysolipid signaling in reproduction, such as ovarian function, fertilization and early embryo development (63). JH III, MF and Met significantly altered membrane phospholipid metabolism, indicating common impacts leading to lysophospholipid G-protein-coupled receptor mediation and activation of related enzymes (64) with a potential activating effect to hepatopancreas. Among all significant metabolites, no juvenoid degradation products were identified. Therefore, we assume the concentration we chose was not excessive, and degradation was not significant.

The relative mRNA expression of genes related to lipid metabolism in response to different hormonal treatments in hepatopancreas of the mud crab were presented in Figure 5 . A significantly upregulated level of genes involved in lipid catabolism was observed in FA treatment as opposed to MF, JH III and Met when compared with control. HNF4 is the evolutionarily conserved transcription factor regulating lipid metabolism to efficiently mobilize stored fat, driving β-oxidation for energy production (65). HNF4 is a direct regulator of β-oxidation enzymes such as CPT (66), congruent with the synchronous expression patterns we found for hnf4 and cpt. In Figures 5A, B , MF, JH III and Met showed similar patterns without significant alterations of hnf4 and cpt expression compared with control, indicative of lipid reserve, whereas FA treatment significantly stimulated lipid oxidation via hnf4 and cpt activation, indicative of lipid utilization. A negative association between HNF4 and TG accumulation, and even oogenesis were reported in adult Drosophila melanogaster (65, 67). In mosquito, JH III decreased HNF4 and lipid oxidation (48). The RNAi depletion of methoprene-tolerant receptor for JH III stimulated TG accumulation in insects (48), demonstrating the possible role of JH III in switching from lipid catabolism to lipid accumulation. MF in crustaceans is similar to JH III in insects. Therefore, the higher TG accumulation in MF, JH III and Met treatments in Figure 4B might result from less lipid oxidated than that in FA treatment. Finally, VLCAD is also involved in fatty acid β-oxidation, but was not affected and did not show any significant difference among all groups in Figure 5C .

In terms of lysophospholipids, the important intermediates in lipid metabolism and ‘lipid signal’ for various other metabolic pathways, we identified an elevated level of LPC 16:0 in FA treatment compared with other treatments and control in Figure 4C . This might be explained by the similarly significantly upregulated ppt expression level in FA group in Figure 5D . PPT is a soluble lipase that cleaves fatty acid chains such as palmitate, and further increases production of LPC 16:0 as catalyzed by cytosolic phospholipase A2. Free fatty acid release from storage lipids may transfer to ovary providing essential nutrients. Binding of saturated fatty acids like palmitoyl CoA (16:0) also results in activation of HNF4 (68) in concert with our data. Moreover, another important role is that PPT lies at the junction of fatty acid elongation and degradation according to KEGG, where hexadecanoyl-CoA are synthesized and then transformed by PPT into hexadecanoate for fatty acid degradation. By initiating lipid degradation during vitellogenesis and oocyte development, which was reported in the brown planthopper (69) and might also be the case of our study, the expression of ppt should subsequently alter the content of total lipids, explaining the differed lipid composition we found in FA group compared with MF, JH III and Met.

On the other hand, KAT is capable of catalyzing the Claisen condensation reaction between the two acyl-CoA, thereby achieving carbon chain elongation (70) and lipid synthesis (71). However, we detected no significant difference among treatments in Figure 5E . Therefore we assumed that the differences in TG accumulation levels were generally due to different levels of lipid degradation. Our results pointed in a direction of lipid retention reflected by the significant accumulated TGs in MF, JH III and Met treatment groups against control, contributing to effective ovarian maturation. However, for FA application, lipid reserves were utilized for energy, and might also be used in nutrient absorption or other signaling pathways. Therefore, JH III and Met shared similar activities with MF, while FA played rather different roles to MF, JH III and Met in mud crab hepatopancreas during ovarian maturation.