Isobaric Tags for Relative and Absolute Quantitation-Based Proteomics Analysis Provides a New Perspective Into Unsynchronized Growth in Kuruma Shrimp (Marsupenaeus japonicus)

Unsynchronized growth is a common phenomenon in farmed crustaceans. The underlying molecular mechanism of unsynchronized growth of crustaceans is unclear. In this study, a comparative proteomic analysis focusing on growth differences was performed using kuruma shrimp Marsupenaeus japonicus, an economic crustacean species, as the model. The study analyzed kuruma shrimp at fast growth stage and steady growth stage from both fast growth group and slow growth group by an Isobaric tags for relative and absolute quantitation (iTRAQ)-based quantitative proteomic analysis method. A total of 1,720 proteins, including 12,291 peptides, were identified. Fifty-two and 70 differentially expressed proteins (DEPs) were identified in the fast growth stage and steady growth stage, respectively. Interestingly, 10 DEPs, including 14-3-3-epsilon-like, GPI, GPD1, MHC-1a, and MHC-1b, were presented in both growth stages. In addition, all these 10 DEPs shared the same expression tendency at these two growth stages. The results indicated that these 10 DEPs are potential growth biomarkers of M. japonicus. Proteins associated with faster growth of M. japonicus may promote cell growth and inhibit cell apoptosis through the Hippo signaling pathway. The fast growth group of M. japonicus may also achieve growth superiority by activating multiple related metabolic pathways, including glycolysis, glycerophospholipid metabolism and Citrate cycle. The present study provides a new perspective to explore the molecular mechanism of unsynchronized growth in crustacean species.


INTRODUCTION
Growth rate is a key factor to associate with the yield and economic benefits for the crustacean farming industry. However, for most crustacean species, the growth is unsynchronized, which means even cultured under the same genetic and environmental conditions (Ren et al., 2017;Jiang et al., 2020;Uengwetwanit et al., 2020), individual crustacean grows at different rates. The unsynchronized growth phenomenon expands the total harvest time, leads to more management efforts and higher farming cost (Hao et al., 2019;Jiang et al., 2020). With more and more large-scale farming models being developed in recent years, how to efficiently minimize the unsynchronized growth becomes a critical question for the industry. However, there is limit information on the underlying mechanism of unsynchronized growth in crustaceans, which hinders the industry development.
Proteomics is one of the most powerful tools available to investigate molecular mechanisms of different phenotypes. Gene expression regulation is a complex multilevel process. In this process, genes are the carriers of genetic information (Cheng et al., 2018), and proteins are the executors of the physiological functions and direct manifestations of life activity (Zhang et al., 2014). However, the correlation between mRNA and protein abundances was only approximately 27∼40% (Maier et al., 2009). Therefore, proteomics analysis focusing on expressed proteins, becomes an effective method to explore the molecular mechanisms corresponding to different phenotypes. Currently, although a number of studies on crustaceans' unsynchronized growth at the gene-level were observed (Lv et al., 2015;Uengwetwanit et al., 2020), there are few related studies focusing on protein identification and quantification (Ren et al., 2017).
iTRAQ-based quantitative proteomic analysis method, a new reliable approach, has been widely used in the field of crustacean proteomics in recent years. iTRAQ allows simultaneous identification and quantification of proteins from multiple samples with high coverage (Ren et al., 2017;Xu et al., 2017). In terms of unsynchronized growth of crustaceans, Ren et al. (2017) identified 30 differentially expressed proteins (DEPs) related to the growth of Portunus trituberculatus, among these 30 DEPs, most of the growth-related proteins were upregulated, including those in metabolism, immune responses, DNA duplication and protein synthesis.
Current studies provide valuable information to better understand the genetic regulation mechanism of unsynchronized growth in crustacean species. However, previous research lacks recognizing the impact of sex and molting cycle, which influence crustacean growth at various levels (Liang et al., 2017;Lin et al., 2019). Unsynchronized growth is a life-long process common for most animals, the study on crustacean will not only illuminate the understanding on invertebrate animals, it will also offer valuable information for our interpretation of the overall biology research.
Kuruma shrimp Marsupenaeus japonicus is a broadly farmed economical shrimp species in the Indo-West Pacific region in both tropical and subtropical waters (Wang et al., 2018(Wang et al., , 2020Liu et al., 2019), where the supply of wild M. japonicus cannot meet the demand. So, shrimp farming industry developed fast in recent years to fill the gap (Jung et al., 2013). In this study, a proteomic analysis comparing the fast growth group and slow growth group of M. japonicus in both fast growth stage (40 days of age with higher specific growth rate) and steady growth stage (70 days of age with lower specific growth rate) by iTRAQ was conducted. The study identified Hippo signaling pathway and some metabolic pathways potentially impact the growth of M. japonicus. Our research not only shed light on molecular mechanisms underlying the growth of M. japonicus but also would promote related research focusing on solve problems in the agricultural production and better advance the farming industry.

Shrimp and Experimental Design
The study protocol was approved by the ethics review board of the Institutional Animal Care and Use Committee of Guangdong Ocean University. Shrimp were obtained from the East Sea Island Marine Research Station of Guangdong Ocean University. The shrimp culture management was as described in our previous report (Zhao et al., 2021). Briefly, the shrimp came from a full-sib family (Family 1) constructed by our team. Ten days after larval metamorphosis to postlarvae (10 days of age), the family was transported to an 800 m 2 ponds. The culture densities was 40 individuals/m 2 . The shrimp were fed 2 times daily with commercial pellet diets with crude protein ≥ 43.0% (Guangdong Yuequn Biotechnology Co., Ltd., China). The zero-exchange water culture pattern was adopted. Throughout the culture experiment, the pond was aerated to maintain the dissolved oxygen at > 5 mg/L. The pH and salinity of the pond remained stable, ranging from 8.16 to 8.62 and 28 to 32 , respectively. The concentrations of ammonia nitrogen, nitrite nitrogen and nitrate nitrogen were kept at low levels, lower than 1.017, 0.206, and 1.336 mg/L, respectively. Based on the growth characterization of kuruma shrimp (Zhao et al., 2021), two time points, 40 days of age (fast growth stage with higher specific growth rate) and 70 days of age (steady growth stage with lower specific growth rate), individuals exhibiting significant differences in size, were selected for subsequent analysis. In addition, since sex is a key factor in determining shrimp development and growth, only female shrimp were selected for the experimental analysis in this study (Lin et al., 2019;Zhao et al., 2021). At each time point, 500 female shrimp were randomly captured. The 10 with the highest body weight were defined as the fast growth group (FG), and the 10 with the lowest body weight were defined as the slow growth group (SG). Then, 3 female shrimp in the intermolt stage (Chan et al., 1988) were sampled in each group. The average body weights of the FG40, SG40, FG70, and SG70 group shrimp were 1.21 ± 0.08, 0.33 ± 0.01, 2.37 ± 0.08, and 0.73 ± 0.03 g, respectively. The body weight of the FG and SG group shrimp at each age was found to be significantly different by t-test. After dissecting each shrimp, the abdominal muscle, the main target organ for shrimp growth, was taken as a sample for analysis (Zhao et al., 2021). There were 3 samples obtained from each group, with 12 samples in total. After dissection, the samples were immediately placed in liquid nitrogen and then stored in an ultralow temperature freezer at −80 • C until protein extraction.

Protein Extraction and Digestion
According to the manufacturer's instructions, 12 muscle samples were extracted using an animal tissue total protein extraction kit (Bangfei Biotech, China). The protein concentration of each sample was quantified using a 2D Quant kit (GE Healthcare, Waukesha, WI, United States). Subsequently, 200 mg of total protein obtained from each sample solution was digested with trypsin (Promega, Madison, WI, United States) at a protein:trypsin ratio of 50:1 at 37 • C for 12 h.

Liquid Chromatography-Mass Spectrometry Analysis
The fractionated peptides were dissolved in buffer A (0.1% FA) and loaded onto a Thermo Scientific reversed-phase C18 precolumn (3 mm, 100 µm × 20 mm). The gradient included 6, 6-28, 28-38, 38-100, and 100% solvent B (0.1% FA and 80% ACN), and the corresponding times were 0-3, 3-45, 45-50, 50-55, and 55-60 min, respectively. The samples were analyzed at a constant flow rate of 300 nl min −1 on an EASY-nLC UPLC system (Thermo Fisher Scientific). The peptides were then analyzed on a Q Exactive TM plus hybrid quadrupole-Orbitrap mass spectrometer (Thermo Fisher Scientific). Intact peptides were detected in the Orbitrap at a resolution of 60,000. Peptides were selected for MS/MS using an NCE setting of 30. Ion fragments were detected in the Orbitrap at a resolution of 15,000. A data-dependent procedure alternating between one MS scan and 20 MS/MS scans was applied for the top 20 precursor ions above a threshold ion count of 10,000 in the MS survey scan with 30.0-s dynamic exclusion. The electrospray voltage applied was 2.0 kV. The m/z scan range was 350-1,800.

Protein Identification and Data Analysis
Protein identification was performed using the Proteome Discover 2.2.0.388 (Thermo Fisher Scientific). The database (SRP305604) used for identification was the M. japonicus transcriptome constructed by our team (Zhao et al., 2021). The parameters for analysis were set as follows: enzyme, trypsin; 15.0 ppm for peptide mass tolerance; up to 2 cleavages; carboxymethyl (C), iTRAQ 8-plex (N-term), iTRAQ 8-plex (K) and iTRAQ 8-plex (Y) as fixed modifications; and oxidation (M) and acetyl (protein N-term) as variable modifications. To reduce the probability of false protein identification, the protein identification was based on at least one unique peptide with 99% confidence and a FDR (false discovery rate) ≤ 0.01. For protein quantification, at least two unique peptides were required. The quantitative protein ratios were weighted and normalized by the median ratio in Mascot. Only proteins with P-values < 0.05 and FC (fold-change) values > 1.0 were considered DEPs. The Gene Ontology (GO) terms 1 and Clusters of Orthologous Groups of Proteins database 2 were used to classify and group all identified proteins. In addition, the COG database and Kyoto Encyclopedia of Genes and Genomes (KEGG) 3 were used to analyze the identified DEPs.

Quantitative Real-Time PCR
To verify the protein expression levels determined by proteomic analysis, the expression levels of myosin heavy chain type-1a (MHC-1a) and myosin heavy chain type-1b (MHC-1b) were further quantified by Quantitative Real-Time PCR (qPCR) using the original protein samples (FG40, SG40, FG70, and SG70). The primers were designed using Primer 5 (Supplementary Table 1), and elongation factor 1-α (Qiao et al., 2017) and 18S ribosomal RNA (Sellars et al., 2007) of M. japonicus were used as internal controls. qPCR was performed in a total volume of 20.0 µl, including 10.0 µl of 2 × SYBR Green PCR Mix, 0.4 µl of each primer (10 µM), 8.2 µl ddH 2 O, and 1.0 µl cDNA (50 ng). The cycling parameters were as follows: 95 • C for 30 s, then 40 cycles of 95 • C for 5 s and holding at 60 • C for 30 s. Dissociation curve analysis was further performed to check if only one PCR product was amplified. Three biological replicates were analyzed for each gene, and three technical replicates were analyzed for each biological replicate. Gene expression levels were compared and converted to fold differences using the comparative CT method (2 − CT ) (Vandesompele et al., 2002).

Protein Profiling
To obtain an overview of the M. japonicus proteome, protein samples were prepared from the fast growth group and slow growth group at the fast growth stage and steady growth stage (FG40, SG40, FG70, and SG70) and quantified by iTRAQ (Figure 1). A total of 1,720 proteins, including 12,291 peptides, were identified. In addition, 1,656 proteins were annotated in at least one of the protein databases, including Nr (1,656), GO (1,301), COG (805), and KEGG (730). In total, 805 proteins were further classified into 25 COG classifications (Figure 2A). The largest category was translation, ribosomal structure, and biogenesis (J) (143), followed by posttranslational modification, protein turnover, chaperones (O) (111), general function prediction only (R) (102), energy production and conversion (C) (68), cytoskeleton (Z) (66), signal transduction mechanism (T) (55), amino acid transport and metabolism (E) (55), and carbohydrate transport and metabolism (G) (51). Moreover, 1,301 proteins were assigned to different classes in three main GO categories, namely, biological process (21 subcategories), cellular component (16 subcategories), and molecular function (12 subcategories) ( Figure 2B). In the category of biological process, the largest number of proteins was assigned to the term cellular process, whereas the binding and cell subcategories contained the largest number of proteins in the categories of molecular function and cellular component, respectively.

Isobaric Tags for Relative and Absolute Quantitation Quantification
In the present study, SG40 vs. FG40 and SG70 vs. FG70 were set as comparison groups. Finally, 52 DEPs (28 upregulated proteins and 24 downregulated proteins) between SG40 and FG40 (Supplementary Table 2  identified ( Figure 3A). Ten DEPs were found to be present in both the SG40 vs. FG40 and SG70 vs. FG70 comparisons ( Figure 3B). Interestingly, all shared DEPs expressed the same tendency in the two comparative groups (Figure 3C). Among these 10 DEPs, glucose-6-phosphate isomerase (GPI) and glycerol-3-phosphate dehydrogenase (GPD1) were upregulated in the fast growth group at both growth stages, whereas the other 8 DEPs, including 14-3-3-epsilon-like and cyclophilin A, were downregulated in the fast growth group at both growth stages.

Differentially Expressed Protein Analysis by COG and Kyoto Encyclopedia of Genes and Genomes
The DEPs in the SG40 vs. FG40 and SG70 vs. FG70 groups were further categorized according to COG classification. Translation, ribosomal structure, and biogenesis, carbohydrate transport and metabolism, and posttranslational modification, protein turnover, and chaperones were the 3 COG classifications enriched with the most proteins in the two comparative groups ( Figure 4A). Moreover, DEP functions and biological processes were analyzed by KEGG. Twenty-five DEPs in SG40 vs. FG40 and 40 DEPs in SG70 vs. FG70 were enriched in 13 and 16 Level 2 pathways, respectively, of which 9 Level 2 pathways were shared by the two comparative groups (Supplementary  Figure 1). Further analysis of the pathways included in these 9 Level 2 pathways revealed 20 pathways in SG40 vs. FG40 and 34 pathways in SG70 vs. FG70, of which 16 pathways were shared by the two comparison groups (Figure 4B). These pathways were involved in the following physiological activities: lipid metabolism (glycerophospholipid metabolism), global and overview maps (carbon metabolism, biosynthesis of amino acids), energy metabolism (oxidative phosphorylation), carbohydrate metabolism (starch sucrose metabolism, pentose phosphate pathway, glycolysis/gluconeogenesis, amino sugar and nucleotide sugar metabolism), amino acid metabolism (tryptophan metabolism, lysine degradation, glycine, serine and threonine metabolism), translation (ribosome), folding, sorting, and degradation (protein processing in the endoplasmic reticulum), signal transduction (the Hippo signaling pathway), and transport and catabolism (phagosome, endocytosis).

Potential Pathways and Genes Related to Unsynchronized Growth
The Hippo signaling pathway is an essential regulator of cell proliferation and apoptosis during development. In this study, we identified some key genes of the Hippo signaling pathway in M. japonicus. As mentioned above, 14-3-3-epsilon-like was downregulated in the fast growth group at both growth stages (Figure 5B), whereas two types of actin1 were upregulated in the fast growth group at the fast growth stage (Figure 5B). In addition, the Hippo signaling pathway had an impact on the downstream pathways. Calcineurin subunit A (CaN) in the Wnt signaling pathway was also upregulated in the fast growth group at the fast growth stage (Figure 5C), and eukaryotic translation initiation factor 5-like (eIF5-like) in the apoptosis pathway was downregulated in the fast growth group at the steady growth stage (Figure 5D). The regulatory network centered on the Hippo signaling pathway, which showed a great impact on muscle growth ( Figure 5A).
Glycolysis plays a central role in the anabolism and catabolism of organisms and provides intermediates for other metabolic pathways and energy for life activities. In this study, both the SG40 vs. FG40 and SG70 vs. FG70 comparison groups had 3 DEPs enriched in this pathway. Among these DEPs, GPI was shared by the two comparison groups ( Figure 6B). Except for phosphoglucomutase (PGM), GPI, 2,3-bisphosphoglycerate-dependent phosphoglycerate mutase (PGAM) and triosephosphate isomerase (TPI) in SG 40 vs. FG 40, and GPI and pyruvate kinase 3 (PK3) in SG 70 vs. FG 70, were upregulated in the fast growth group (Figure 6B). In addition, a variety of metabolic pathways related to glycolysis were affected. Fructose 1,6-bisphosphatase (FBP) in the gluconeogenesis (Figure 6C), GPD1 in the glycerophospholipid metabolism pathway ( Figure 6D) and 2-oxoglutarate dehydrogenase (OGDH) in the citrate cycle were all upregulated in the fast growth group (Figure 6E). These pathways and glycolysis form a complex energy conversion network ( Figure 6A).

Quantitative Real-Time PCR Validation
MHC-1a and MHC-1b were verified by qPCR, and the qPCR results were consistent with the proteomic analysis results (Figure 7). MHCs are related to the growth of crustaceans (Jung et al., 2013). Compared with other MHCs identified in this study, the expression levels of MHC-1a and MHC-1b in the fast growth group were lower than those in the slow growth group during both the fast growth stage and steady growth stage.

DISCUSSION
In this study, a comparative proteomic analysis of M. japonicus between fast growth group and slow growth group in both fast growth stage and steady growth stage by iTRAQ was conducted. The study identified some proteins related to growth, Hippo signaling pathway or metabolic pathways potentially impact unsynchronized growth of M. japonicus, which is the first report on identification and quantification of proteins on this topic with the iTRAQ method.
The proteins involved in unsynchronized growth mechanism of M. japonicus in both fast growth stage and steady growth stages were consistent. Sixteen pathways related to lipid metabolism, global and overview maps, energy metabolism, carbohydrate metabolism, amino acid metabolism, translation, folding, sorting and degradation, signal transduction, and transport and catabolism, were shared by the two different growth stages. The data revealed that individual growth differences were the results of the continuously development of these physiological activities. Among these proteins, 10 DEPs showed consistent expression tendency in both fast and slow growth stages, which are more likely to weigh in and regulate the growth of M. japonicus. Further analysis concentrating on these 10 DEPs led us to identify a Hippo signaling pathway.
Hippo signaling pathway is a signaling pathway that controls organ size in animals through the regulation of cell proliferation and apoptosis (Tapon et al., 2002;Harvey et al., 2003;Jia et al., 2003;Pantalacci et al., 2003;Thompson and Cohen, 2006). A three-step kinase cascade trigger the center reactions in the Hippo signaling pathway and a transcriptional co-activator, Yorkie (Yki) is the key effector for the kinase cascade. Yki can interact with transcription factors in the nucleus and lead to the transcriptional upregulation of a number of target genes (Wang et al., 2016;Liang et al., 2018). Cytoskeletal protein 14-3-3 is another protein that can interact with phosphorylated Yki and cause Yki being detained in the cytoplasm (Dong et al., 2007;Oh and Irvine, 2008). When 14-3-3 is inactivated, it enhances tissue overgrowth which was induced by Yki's overexpression (Ren et al., 2010). Remarkably, a similar regulation was observed in our study. Shrimps with a low expression level of 14-3-3 (14-3-3epsilon-like) presented growth superiority. The results confirmed that the Hippo signaling pathway influenced the growth in M. japonicus. The Hippo signaling pathway had an impact on downstream pathway genes. For example, gene CaN, which is related to cell proliferation (Chen et al., 2003) was upregulated in the fast growth group, and the gene eIF5-like related to cell apoptosis (Li et al., 2004) was downregulated. The data signals that the fast growth group of M. japonicus achieved a growth superiority by promoting cell growth and inhibiting cell apoptosis through the Hippo signaling pathway.
Another notable finding is that almost all the differentially expressed glycolytic proteins (GPI, PGAM, TPI, PGM, PK3) were upregulated in the muscle of fast growth group of M. japonicus, which indicates that glycolysis assists in the growth superiority of M. japonicus. Glycolysis is related to muscle development. Mutation of glycolytic genes can lead to different muscle disease (García et al., 2009;Tixier et al., 2013). In return, muscle diseases can also result in glycolytic enzyme activity decrease and glycogen accumulation (Sharma et al., 2003;Coley et al., 2012). In fishes, almost all glycolysis genes in muscles were found to be upregulated in hybrid groupers showing growth superiority (Sun et al., 2016).
In addition, FBP in gluconeogenesis, GPD1 in glycerophospholipid metabolism and OGDH in Citrate cycle may also promote growth superiority in M. japonicus, since all these proteins were upregulated in the fast growth group. The reason behind the activated metabolic pathways pointed to the interaction between the glycolysis and metabolic pathways: Certain metabolites in glycolysis also are intermediates in other anabolic pathways. Therefore, the activation of the glycolysis led to the activation of multiple metabolic pathways, including gluconeogenesis (Cohen, 2014), glycerophospholipid metabolism and Citrate cycle (Panchal et al., 2001;Sharma et al., 2005).
Furthermore, MHC-1a and MHC-1b inhibits growth superiority in M. japonicus, as they were downregulated in the fast growth group, confirmed by both proteomics and qPCR methods. Our data verified a previous study on a 60-day-old Siniperca chuatsi . MHC contributes to muscle growth. There are studies showing that the expression of MHC correlates with the growth rate of vertebrates (Dhillon et al., 2008;Churova et al., 2015;Overturf and Hardy, 2015). In crustaceans, the expression of myosin and actin correlated with both abdominal and skeletal muscle morphology changes and biochemical composition changes in different molt stages (Medler, 2003;Cesar et al., 2006). During myogenesis, MHC binds firmly to muscle precursor cells to promote muscle development (Kreissl et al., 2008;Harzsch and Kreissl, 2010;Jirikowski et al., 2010). MHC is a large family including many members with different functions, so how to further characterize MHC-1a and MHC-1b's functions related to the growth of M. japonicus will be an excited topic for future studies.

CONCLUSION
In conclusion, the present study for the first time investigated the mechanism of unsynchronized growth in shrimp by iTRAQ proteomics analysis. In our study, we analyzed DEPs between the fast growth group and slow growth group in different growth stages in M. japonicus and identified 10 key DEPs, including 14-3-3-epsilon-like, GPI, GPD1, MHC-1a and MHC-1b could impact M. japonicus's growth. The shrimp M. japonicus can achieve a growth superiority by activating one or more pathways, such as Hippo signaling pathway and metabolic pathways with glycolysis. Our study provides a new perspective for understanding the mechanism of unsynchronized growth of crustacean species.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are publicly available. This data can be found here: ProteomeXchange Consortium (http://proteomecentral.proteomexchange.org) via the iProX partner repository with the dataset identifier PXD028033.

AUTHOR CONTRIBUTIONS
JZ, CS, and ZN conceived and designed the study. JZ performed the experiments, analyzed the data, and wrote the manuscript. ML, ZL, YH, YZ, LL, and GC assisted the experiments. CS and ZN revised the manuscript. All authors contributed to the article and approved the submitted version.