Whole Transcriptome Analysis Provides Insights into Molecular Mechanisms for Molting in Litopenaeus vannamei

Molting is one of the most important biological processes in shrimp growth and development. All shrimp undergo cyclic molting periodically to shed and replace their exoskeletons. This process is essential for growth, metamorphosis, and reproduction in shrimp. However, the molecular mechanisms underlying shrimp molting remain poorly understood. In this study, we investigated global expression changes in the transcriptomes of the Pacific white shrimp, Litopenaeus vannamei, the most commonly cultured shrimp species worldwide. The transcriptome of whole L. vannamei was investigated by RNA-sequencing (RNA-seq) throughout the molting cycle, including the inter-molt (C), pre-molt (D0, D1, D2, D3, D4), and post-molt (P1 and P2) stages, and 93,756 unigenes were identified. Among these genes, we identified 5,117 genes differentially expressed (log2ratio ≥1 and FDR ≤0.001) in adjacent molt stages. The results were compared against the National Center for Biotechnology Information (NCBI) non-redundant protein/nucleotide sequence database, Swiss-Prot, PFAM database, the Gene Ontology database, and the Kyoto Encyclopedia of Genes and Genomes database in order to annotate gene descriptions, associate them with gene ontology terms, and assign them to pathways. The expression patterns for genes involved in several molecular events critical for molting, such as hormone regulation, triggering events, implementation phases, skelemin, immune responses were characterized and considered as mechanisms underlying molting in L. vannamei. Comparisons with transcriptomic analyses in other arthropods were also performed. The characterization of major transcriptional changes in genes involved in the molting cycle provides candidates for future investigation of the molecular mechanisms. The data generated in this study will serve as an important transcriptomic resource for the shrimp research community to facilitate gene and genome annotation and to characterize key molecular processes underlying shrimp development.

Introduction findings, our understanding of the molecular mechanisms underlying shrimp molting remains very limited.
In order to connect the mechanisms and molecular events associated with shrimp molting, we used RNA-Sequencing (RNA-seq) to investigate expression changes across all genes during the molting process of L. vannamei, one of the most commonly cultured shrimp species, representing nearly 40% of penaeid shrimp production worldwide [27]. This approach is considered to have high power to dissect gene networks associated with particular biological and developmental processes at a whole transcriptome level [28][29][30]. Here, we report a comprehensive analysis of global transcript expression changes associated with shrimp molting, which not only enhances our understanding of the molecular mechanisms underlying shrimp molting, but also provides a valuable resource of transcriptome data shrimp genome annotation and further identification of candidate genes controlling molting in shrimp and other molting animals.

Sample collection and RNA isolation
Healthy adult Pacific white shrimp (L. vannamei) with an average body length of 14-16 cm were collected from laboratory culture ponds. All shrimp were the same generation and had been cultured over six months to minimize generational and environmental effects. All procedures involving animals throughout the experiments were conducted in strict accordance with the Chinese Legislation on the Use and Care of Laboratory Animals. And all animal experiments were performed as per the institutional ethic committee guidelines of Institute of Oceanology Chinese Academy of Sciences (IOCAS).
Animals were sorted by molt stage according to the appearance of epidermis, pigmentation, formation of new setae, and presence of matrix or internal cones in the setal lumen (Fig 1). Molting stages were classified as inter-molt (C), pre-molt (D0, D1, D2, D3, and D4), and postmolt (P1 and P2). A loop design comparing consecutive molt stages (C-D0-D1-D2-D3-D4-P1-P2-C) was employed for the analysis to achieve consecutive times for all periods without interruption. Each whole shrimp was considered a separate sample and two samples were taken as biological replicates for each molt stage. All samples were immediately frozen in liquid nitrogen and stored at -80°C until RNA isolation.
Each sample was crushed in a mortar with liquid nitrogen, and total RNA was prepared using a Trans-up RNA isolation kit (Biostar, Shanghai, China) according to the manufacturer's protocol. The yield and purity of each RNA sample were determined using a NanoDrop™ 2000 spectrophotometer (Thermo Scientific, USA), and the integrity of all RNA samples was assessed by gel electrophoresis with 1% agarose. Total RNA was treated with DNase to remove DNA contamination.

RNA sequencing library construction and Illumina sequencing
Isolation and enrichment of mRNA from total RNA was performed using oligo (dT) magnetic beads (Illumina, CA, USA). Then, mRNA was fragmented to short fragments to be used as templates for random hexamer-primed synthesis of first strand cDNA by fragmentation buffer. Second-strand cDNA was synthesized using buffer, dNTPs, RNase H, and DNA polymerase I. A paired-end cDNA library was synthesized using the Genomic Sample Preparation Kit (Illumina, CA, USA) according to the manufacturer's instructions. Short fragments were purified with QIA-Quick 1 PCR extraction kit (QIAGEN, Germany) and eluted in 10 μL of EB buffer (QIAGEN, Germany). These short fragments were connected via sequencing adapters (Illumina, CA, USA). Agarose gel electrophoresis was used to select fragments approximately 50 bp in size. Finally, cDNA libraries were sequenced on an Illumina HiSeq™ 2500 (Novogene, Beijing).

Assembly of sequencing data and gene annotation
Raw sequence data were transformed by base calling into sequence data and stored in fastq format. Raw reads were cleaned by removing adapter sequences, empty reads, and low quality sequences. Cleaned reads were assembled with the software package RSEM [31] using L. vannamei reference transcriptome data obtained from the molting-transcriptome sequencing (SRX1411196) and combined with the data previously sequenced by our laboratory (SRR1460493, SRR1460494, SRR1460495, SRR1460504 and SRR1460505).
For annotation analysis, unigenes were BLASTX-searched against five databases, including the National Center for Biotechnology Information (NCBI) non-redundant protein sequence (NR) database, the NCBI non-redundant nucleotide sequence (NT) database, KEGG Orthology (KO) database, Swissprot, and the PFAM database, using a cut-off E-value of 10 −5 . Unigenes were annotated based on BLASTX results, and the best alignments were used for downstream analyses.

Normalized expression levels of genes from RNA-seq
To eliminate the influence of different gene lengths and sequence discrepancies on expression calculations, gene expression levels based on read counts obtained by RSEM were normalized using the FPKM (Fragments Per Kilo bases per Million fragments) transformation [32]. Thus, calculated gene expression levels could be used for direct comparison among samples. Expression values were standardized across the dataset to enable the data from different genes to be combined.

Screening of differential expression genes (DEGs)
Using the R package DEGseq, differentially expression genes (DEGs) were identified with a random sampling model based on the read count for each gene at different developmental stages [33]. False discovery rate (FDR) 0.001 and absolute value of log 2 ratio 1 were set as the threshold for significance of gene expression differences between adjacent samples (C-D0, D0-D1, D1-D2, D2-D3, D3-D4, D4-P1, P1-P2, and P2-C).

Gene Ontology and KEGG analysis
Gene Ontology (GO) terms were used to describe biological processes, molecular functions, and cellular components. As an international standardized gene functional classification system, GO offers both a dynamically updated controlled vocabulary and strictly defined concepts to describe the properties of genes and their products comprehensively. The Blast2GO(version 3.0) (https://www.blast2go.com/) program was used to obtain GO annotations for all genes with a Fisher's Exact Test (filtered with FDR 0.01) [34]. The unigene sequences were aligned to the Clusters of Orthologous Group (COG) database. Using GO functional classification analysis (WEGO), we categorized all genes based on function [35]. The Kyoto Encyclopedia of Genes and Genomes (KEGG) database was used to assign and predict putative functions and pathways associated with the assembled sequences [36]. A heat map which grouped genes according to FPKM values was generated in Cluster 3.0 [37] and visualized in TreeView 1.6 to analyze expression levels across molting periods [38].

Real time qPCR amplification
To validate RNA-seq data and expression profiles, six genes were randomly selected for validation using real-time quantitative polymerase chain reaction (RT-qPCR). Actin T2 (c82047_g1) was used as an internal standard, and relative gene expression levels were calculated using the comparative Ct method with the formula 2 -ΔΔCt [39]. All samples were run in triplicate in separate tubes; each cDNA sample was run in duplicate. All data were expressed as mean +SD after normalization. Real-time qPCR results were then compared with transcriptome data to detect the expression correlation of each gene. The primers used for amplification and the annotations of the products are listed in S1 Table. Data Availability The sequence data in this study have been deposited into the NCBI Sequence Read Archive (http://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?view=studies), and the accession numbers of the eight SRA samples are as follows: SRX1098368, SRX1098369, SRX1098370, SRX1098371, SRX1098372, SRX1098373, SRX1098374, and SRX1098375.

Result of RNA-Sequencing and analysis
Sixteen normalized cDNA libraries prepared from all molting stages were sequenced using the Illumina Hiseq 2500 platform. After removal of adaptor sequences, duplicate sequences, ambiguous reads, and low-quality reads, a total of 148.25 millions of clean reads (7.41Gb) with 98.7% Q20 bases were selected as high-quality reads for expression analysis. An overview of the sequencing results is summarized in Table 1. Among the total number of reads from sixteen samples, 83.11% to 88.38% were matched in comparison with the reference genome. The percentages of mapped reads were similar in these libraries.
By comparing with the reference transcriptome data, all clean reads were assembled into 93,756 unigenes. To understand potential genetic mechanisms undergirding L. vannamei molting, 93,756 assembled unigenes were BLASTX-searched against five databases. Of the 93,756 unigenes, 15,582 match to known proteins in the NR database, and 3,689 match to putative homologues in the NT database. The KO database provided annotation for 6,493 unigenes, the Swissprot database confirmed matches for 12,873 unigenes, and 20,784 unigenes found putative homologues in the PFAM database.  , respectively, were found to be expressed stage-specifically. In contrast, a large number (28,508) of unigenes were found to be expressed at all eight stages. Some genes exhibited low variance or low expression levels, and were thought either to be housekeeping genes or expressed infrequently during the molting cycle.

Analysis of differentially expressed genes (DEGs) among molting stages
Our analysis targeted genes expressed at relatively high levels between adjacent molting stages. Using log 2 ratio 1 and FDR 0.001 as a threshold, a total of 5,117 differentially expressed genes (DEGs) between any two adjacent molting stages (i.e., C-D0, D0-D1, D1-D2, D2-D3, D3-D4, D4-P1, P1-P2, and P2-C) were identified (S3 Table). The numbers of up-regulated and down-regulated DEGs are shown in Fig 2. A total of 4,015 and 1,456 DEGs were detected between D4 and P1 and between P1 and P2, respectively, representing the two largest groups of DEGs among the eight comparisons.

Hierarchical clustering of DEGs among the eight stages
Using log ratio values, we performed hierarchical clustering of DEGs expression. Expression levels during the molting cycle were divided into twenty-four categories based on K-means clustering. Detailed expression profile clusters during the eight stages are shown in S1 Fig (Fig 3). These expression patterns not only indicate the diverse and complex interactions among genes, but also suggest that unigenes with similar expression patterns may have similar functions in the molting cycle.

Gene ontology (GO) analysis of RNA-seq data
To associate genes exhibiting different expression patterns with morphological and physiological changes during the molting cycle, we performed gene ontology (GO) enrichment analysis using BLAST2GO. A total of 11,501 GO terms were associated with all DEGs. According to the secondary classification of the GO terms, these DEGs were sorted into 47 functional groups that belong to three main GO categories: biological processes, cellular components, and molecular functions, comprised of 4,885 (42.5%), 3,606 (31.4%), and 3,010 (26.1%) DEGs, respectively (Fig 4). In the three categories, the six subcategories; cell (GO:0005623), cell part (GO:0044464), binding (GO:0005488), catalytic (GO:0003824), cellular process (GO:0009987)

Expression pattern of hormone regulation genes
Hormones have an important role in controlling the crustacean molting cycle. In our study, eleven genes related to MIH and crustacean hyperglycemic hormone (CHH) were identified. Interestingly, except for two unexpressed genes, down-regulated expressions were detected in D2-D3 in all others (Table 3). Moreover, we summarized the variance in transcript levels of molting related hormones. As a result, thirteen types of molting factors related to hormone regulation were identified. Interestingly, a trend of up-regulated transcript levels was found in all thirteen factors during D3-D4 (Table 4).

Expression profile of skelemin genes
Expression analysis revealed high transcript levels of myosin, troponin, and actin following molting behavior (D4-P1). Members with high expression (FPKM value >1000 in at least one stage) were grouped together, and 49 members were found (6 were actin, 39 were myosin and 4 were troponin) ( Table 5). Interestingly, 98% of the factors showed an up-regulated expression in D4-P1 (Table 5 bolded text). Some factors showed large changes in expression. For example, the FPKM values of actin T2 (c82047_g1) and skeletal muscle actin 6 (c78558_g1) increased more than 11,000 at D4-P1. Additionally, the FPKM values of myosin light chain (c26953_g1) and myosin light chain 2 (c61912_g2) increased more than 7,800, and those of troponin T (c70761_g2 and c59125_g1) increased more than 3,300 during molting.

Expression pattern of immune related genes
To reveal the expression changes of the immune-related genes during molting, we summed the expression patterns of the immune pathways as a proxy for the immunity level of shrimp according to Li and Xiang (2013a,b) [53,54]. Thirty-six immunity families, including 236 immune genes, were categorized between molting stages. Clusters displayed expression Expression changes of hemocyanin, chitinase and serine protease superfamily in molt process To identify the possible factors involved in the molting process, we investigated the expression profile of members of the hemocyanin, chitinase and serine protease superfamily during the molting process. Data mining of the annotated genes identified twelve genes associated with hemocyanin. Similarly, 19, 16, and five genes were found to be related to serine proteinase, trypsin, and chymotrypsin, respectively. Twenty-three unigenes related to chitinase were  (14) Lysosome (9024) Regulation of actin cytoskeleton (-23032) 6 Peroxisome (17) Cardiac muscle contraction (16315) Antigen processing and presentation (-7150) Starch and sucrose metabolism (13) Ribosome (7665) Cardiac muscle contraction (-22476) 7 Other glycan degradation (16) Phototransduction-fly identified in this study. The expression profile of hemocyanin genes showed a characteristic pattern of up-regulation at inter-molt and down-regulation during D1-D3. In D3-D4, expression was up-regulated, and then declined to the minimum in P1 (Fig 6). Surprisingly, the expression profile of hemocyanin and chitinase were consistent. The expression data revealed that genes corresponding to chitinase were preferentially expressed during periods D1 and D4. Moreover, similar expression patterns were observed in three serine protease superfamily members, including serine protease, trypsin, and chymotrypsin. Serine protease, trypsin, and chymotrypsin displayed expression changes similar to hemocyanin and chitinase, up-regulated in C-D0, D3-D4, and P1-P2 (Fig 6).

Real-time quantitative PCR validation of RNA-seq results
To validate sequencing data, six differentially expressed genes were selected for real-time qPCR analysis. Expression patterns of the selected genes determined by real-time qPCR (Fig 7) are consistent with those determined by RNA-seq, corroborating our results.

Establishment of an overall transcriptome database for the molting process
The results of the transcriptome analysis presented here establish a genome-wide transcriptional landscape for further studies on the developmental and molecular aspects of the molting process in L. vannamei. In terms of RNA-seq results, it was found that the expression changed markedly in D4-P1 and P1-P2. It is known that molting behavior occurs in D4-P1, while P1-P2 is the recovery stage [2]. There is an obvious morphological diversity between these stages. However, the molecular mechanisms involved in molting behavior are not fully elucidated. Our RNA-seq results provide useful insights into the molecular processes underlying these behaviors. By profiling gene expressions and expression changes, it was found that down-regulated genes were predominant in D4-P1, while in P1-P2, various genes were mainly up-regulated  (Fig 2). Based on this result, the DEGs enriched in GO terms related to the shrimp metabolism level were significantly down-regulated during D4-P1. The top five down-regulated terms were GO:0009987 "cellular process" (294 down-regulated DEGs), GO:0008152 "metabolic process" (288 down-regulated DEGs), GO:0005488 "binding" (276 down-regulated DEGs), GO:0071704 "organic substance metabolic process" (249 down-regulated DEGs), and GO:0003824 "catalytic activity" (232 down-regulated DEGs). While the GO terms in P1-P2 were majorly up-regulated, the terms which showed down-regulation in D4-P1 were up-regulated back. Moreover, DEGs for pancreatic secretion (map04972), protein digestion and absorption (map04974), neuroactive ligand-receptor interaction (map04080), lysosome (map04142), and other pathways were down-regulated during D4-P1. In P1-P2, however,  protein digestion and absorption, neuroactive ligand-receptor interaction, and lysosomerelated pathways were recovered (up-regulated). These results were consistent with the morphological changes and the variation in molecular changes between molting behavior and post-molt. In summary, the RNA-seq results provide the basal understanding for the molecular changes occurring during the molting process. Hormone regulation process in L. vannamei molting It is now well established that MIH and CHH peptide families serve as key regulators of hormones controlling the molting process [55]. MIH, which is secreted from a neuro-secretory  center termed the X-organ/sinus gland complex (XO/SG) in crustacean eyestalks, is a member of the CHH family [13,56]. Extirpation of XO/SG can reduce MIH secretion, enable ecdysteroid synthesis, and improve molt frequency [57]. In our study, down-regulated expressions of all MIH and CHH were detected in D2-D3 (Table 3). Based on the inhibitory effect of MIH and CHH on ecdysone secretion, this result suggests that down-regulation of MIH and CHHrelated genes may trigger the molting process, and D2-D3 likely represents the onset of molt hormone-regulation. This speculation is supported further by physiological evidence that reduction in sinus gland MIH content occurs during the late pre-molt stage in Procambarus clarkia [58].The concentration of MIH was found to decrease in pre-molt, followed by an increase in ecdysone in the hemolymph [55]. Furthermore, qualitative measurements of MIH mRNA transcript abundance found MIH peptide synthesis to be dramatically reduced during late pre-molt in eyestalks of Callinectes sapidus [59].
Only few studies have focused on the molecular role of molting-related hormones followed by MIH and CHH. To investigate hormones downstream of MIH and CHH in the molting process pathway and elucidate possible roles in this process, we summarized the variance in transcript levels of molting related hormones. According to the RNA-seq results, a trend of up-regulated transcript levels was found in all thirteen factors during D3-D4 (Table 4). Most of these factors were reported previously to be associated with molting; for example, ecdysone receptors (ECRs) and retinoid-X receptors (RXRs) have been found to be potential targets for hormonal control during limb regeneration and molting in Uca pugilator [45]. In this study, we also identified two ECRs genes (c72141_g9 and c75175_g1) and two RXRs genes (c74810_g1 and c74810_g2). Methyl farnesoate (MF) and farnesoic acid (FA), secreted by the mandibular organs (MO), might have specific stimulatory effects on Y-organs where ecdysteroids are secreted [47]. Synthesis of the vitellogenin gene (VG) in Cherax quadricarinatus has been shown to be induced by 20E, it has been suggested that VG is an ecdysteroid-responsive gene in the molt process [49]. In this study, a trend of up-regulated transcript levels was found in all twelve factors during D3-D4 (Table 4). In light of the fact that the D3-D4 stage occurs immediately before ecdysis, this result suggests that these hormones are associated with molting during the D3-D4 stage. Indeed, D3-D4 likely serves as the implementation phase for molt hormone-regulation. Moreover, the presence of specific up-regulated genes in the molt phase indicates that they could participate in hormone regulation and promote the molting process, as seen in three ecdysteroid-regulated genes (c19570_g1, c81580_g1, and c65832_g1) whose expression increased 2.31, 2.82, and 3.45 fold during D3 to D4, respectively. Similarly, two molting fluid carboxy peptidase precursors (c80306_g1 and c74916_g2) increased 3.21 and 24.5 fold in D3-D4. These genes, potentially involved in hormone regulation, will serve as promising candidates for future analyses.
In summary, we infer that hormone-regulated molting signals likely play distinctly different roles in shrimp molting. Our findings show the first step of molting to be the down-regulation of MIH and CHH family peptides in D2-D3 to initiate the process. Once inhibition by MIH and CHH is prevented, downstream hormones and other factors can participate in regulation in the D3-D4 stage of molting.

Analysis of skelemin in molting process
Post-molt is a period of rapid enlargement and growth for shrimp as a result of body size expansion and exoskeleton reconstruction after ecdysis. The process is complex, involving not just fast water absorption, but also expansion of the skelemin to provide a scaffold and muscle to fill the new body [60]. Skelemin is important for body reconstruction after molting [61]. However, factors participating in enlargement as skelemin during the molting process have remained largely unknown. According to the RNA-seq results, genes related to actin, myosin and troponin were markedly up-regulated after molting (D4-P1) ( Table 5). Actin is a globular multi-functional protein that forms microfilaments and serves as the major component of muscle [62]. It is reported that the transcriptional response of actin (PotActinSK1, PotActinSK2, PotActinHT, and PotAc-tinCT) to 20-hydroxyecdysone (20E) injection was different in Portunus trituberculatus, suggesting that 20E affects muscle plasticity [63]. Myosin is the major component of thick filaments in muscle, and troponin is a complex of three regulatory proteins (troponin C, troponin I, and troponin T) integral to muscle contraction in skeletal and cardiac muscle [64]. As the major contractile protein in vertebrates and invertebrates, myosin was identified to be associated with muscle atrophy during molting, which was characterized by decrease in fiber width and myofibril cross-sectional areas, increase in inter fibrillar spaces, and degradation [65,66]. Troponin and tropomyosin are modulatory proteins that modulate the binding of actin to myosin [67]. Indeed, body enlargement after molting requires frames and skelemin. Considering the expression profiles of actin, myosin, and troponin after molting (Table 5), our findings suggest that these three proteins are expressed abundantly along with skelemin to support the new body of the shrimp after the phase of molting behavior is complete.
Immunity tactic during molting cycle in L. vannamei Immediately following ecdysis, the new exoskeleton is soft and the exercise capacity of shrimp is very weak, so post-molt (P1) is the most sensitive period to suffer an attack from bacteria, viruses, or other predation. It is known that shrimp practice cannibalism which can abet the outbreak of infectious diseases such as White Spot Syndrome (WSS) [68] and the targets of cannibalism are mainly post-molt shrimp. The immune and protection tactics that allow post-molt shrimp to avoid stress are unknown. To address questions regarding post-molt defense, we summed the expression patterns of the immune pathways as a proxy for the immunity level of shrimp according to Li and Xiang (2013a,b) [53,54]. Interestingly, most immune factors were up-regulated in D3-D4. And during molting behavior (D4-P1), parts of immune factors were up-regulated, and the majority up-regulated factors were concentrated in the some downstream and executive factors of immune pathway, like crustin and ALF. Invertebrates including shrimp have no real adaptive immunity and mainly depend on their innate immunity against a variety of pathogens [69]. Antimicrobial peptides (AMPs) are important effectors in innate immunity, which is the front line of host defense against infection by microbes including bacteria, fungi, and viruses [70,71]. Crustin and ALF are important AMPs that play an important role in the immunity of shrimp [72]. Crustin (MjCru I-1) found in Marsupenaeus japonicuscan bind to the cell wall molecules of bacteria, such as lipopolysaccharide (LPS), peptidoglycan (PGN), and lipoteichoic acid (LTA) [69]. LvALF from Litopenaeus vannamei, MjALF from Marsupenaeus japonicus, PtALF from Portunus trituberculatus, and EsALF from Eriocheir sinensis have also been reported to show broadspectrum activity against gram-positive and gram-negative bacteria, fungi, and viruses [73][74][75][76].
In summary, we propose that shrimp will up-regulate the expression of immune-related genes initially in D3-D4. After ecdysis, several downstream factors, such as the antimicrobial peptides, and other factors, such as the i-type lysozyme-like protein 2 (c72940_g1) that is upregulated about 25-folds between D4 and P1, are likely to be expressed to generate an immune response. Once the exoskeleton hardens, the shrimp acquire new protection, and the expression of immune genes recovers to the baseline level.
Expression changes of hemocyanin, chitinase and serine protease superfamily in molt process Some factors have been reported to be associated with the molting process in previous studies. It reported that hemocyanin highly expressed in both the inter-molt and pre-molt periods is reflective of the dual functionality of hemocyanin in preparation for arthropod ecdysis in Portunus pelagicus [77]. Studies have also confirmed that the molting process can affect hemocyanin levels in L. vannamei juveniles [78]. Chitinase is the main enzyme hydrolyzing the glycosidic bonds in chitin to digest old exoskeleton partially [79]. The structures of serine protease superfamily are very similar, even though they recognize different substrates [80][81][82]. The serine protease superfamily members and hemocyanin share phenoloxidase (PO) activity, and studies on L. vannamei have suggested the regulation of trypsin biosynthesis may be, at least in part, under the influence of ecdysteroid hormones [83]. Hence, in our study, expression of serine protease superfamily members correlates strongly with chitinase expression. Taken together with the functional relatedness of these enzymes, this result suggests a shared or similar function in degradation or digestion in the molting cycle. Moreover, the similarity in gene expression profiles of hemocyanin and these three members of the serine protease superfamily may be indicative of activation of the PO pathway. In this study, the RNA-seq results also illustrated that these factors might be associated with molting in L. vannamei.
Other factors that may be associated with the molting process Aquaporins (AQPs) are essential to facilitate the transport of water and other small polar molecules across cell membranes and play an important role in osmotic regulation [84]. After their escape from the confines of a cuticle, shrimp need to rapidly take up water [2]. One aquaporin gene, c76592_g1, was identified in the molting transcriptome and its expression level was the highest after molting behavior (P1), which indicated that it played a role in osmotic regulation during the molting process.
The formation of chitin is catalyzed by chitin synthetase, a highly conserved enzyme involved in chitin synthesis [79]. One chitin synthetase gene (c79847_g1) was found in the transcriptome result; expression analysis showed that it reached a peak in D4 and P2. Along with its functional annotation, chitin synthetase has been suggested to be involved in new cuticle synthesis in D4 and P2.
Myostatin (MSTN/GDF11) was identified to act as a negative regulator of muscle development in vertebrates and in the shrimp Macrobrachium nipponense. During the molt cycle, the expression of Mn-MSTN/GDF11 mRNA was up-regulated significantly at the early post-molt stage, but later decreased gradually [85]. We also found one myostatin gene (c72168_g1) in the molting transcriptome, and the expression pattern was similar to that in Macrobrachium nipponense, with an up-regulation at post-molt, followed by a down-regulation. These results indicated that myostatin may play a role in the molt cycle.
An earlier report showed that knockdown of a novel G-protein pathway suppressor 2 (GPS2) (GenBank accession number JN714124) led to mortality by exuvial entrapment during ecdysis in shrimp (Penaeus vannamei) [86]. In our study, we identified a GPS2 gene, c74297_g2, and its expression exhibited low variance levels; this gene was thought to be a housekeeping gene that acts during the molting cycle in L. vannamei.

Comparison with results of transcriptomic analyses in other arthropods
There have been a number of previous reports on the transcriptomic analysis of the molting process in other arthropods. By analyzing the gene expression profiles in whole animals and organs previously identified to be related to molting by microarray analysis, Kuballa and Elizur performed differential expression profiling of components associated with exoskeletal hardening in the crab P. pelagicus. Moreover, it was found that the C-type lectin receptor, mannose binding protein, PO (trypsin-like, clotting protein precursor-like and antimicrobial proteins) inducers and antimicrobial proteins, in conjunction with hemocyanin, displayed molt cycle-related differential expression profiles, which indicated that they had possible regulatory functions in the calcification or sclerotization cascade of the cuticle [87]. A microarray-based investigation was performed using mRNA collected from the larval stage to the pupal stage for analyzing the cuticular genes in the silkworm, Bombyx mori [88]. A total of 227 cuticular protein genes were found to be expressed, and their expression profiles and clustering was discussed. Custom cDNA microarrays were constructed for the crab Portunus pelagicus to analyze specific differential genes in the molt cycle [77]. By analyzing the expression profiles in thoracic exoskeleton formation during the pupal-to-adult molting in Apis mellifera by microarray analysis, it was found that 1,253 unique DEGs and 547 were up-regulated in the thoracic dorsum after molting, which suggested an induction of expression by the ecdysteroid pulse [89]. Furthermore, 28 of the 36 musclerelated DEGs were up-regulated during the formation of striated fibers attached to the exoskeleton. Another transcriptomic analysis focused on chitin metabolism in four stages (inter-molt, early pre-molt, late pre-molt and post-molt) in the crayfish Cherax quadricarinatus with regard to three key chitin metabolic processes, chitin synthesis, chitin breakdown, and the junction between the metabolism of simple sugars and amino sugars [90].
There were some similarities compared to the results of our work, for example, some common factors are considered to be related to the molt cycle, like the cuticular genes, PO family, chitin synthase, antimicrobial proteins and hemocyanin. However, several differences exist between our results and those of previous studies. These differences may be due to the difference in the techniques used (RNA-seq and microarray). The transcriptome database obtained using Illumina RNA-seq in our study was larger than that obtained by microarray-based analysis in other molting animals, and more genes were detected at each stage of the molting process in our study. Secondly, previous studies mainly focused on exoskeleton reconstruction, but hormone regulation, expression profiles during molting behavior, skelemin, and immunity tactics have not been analyzed so far. We focused on these issues and provided insights into the molecular mechanisms underlying the molting process. This is also the first time that the most important economic species of shrimp, L. vannamei, was used to analyze the transcriptome of whole animals during the molting cycle.
We obtained a considerable number of expressed genes; however, we did not investigate the expression of all genes and focused only on the highly expressed ones. Therefore, some infrequently expressed factors that play important roles in molting may be absent. This limitation can be solved by analyzing the transcriptomes of different tissues and developmental stages in shrimp [91][92][93][94][95] or genes with limited expression.

Conclusion
Here, we report the first comprehensive study utilizing RNA-seq to characterize genes associated with shrimp molting and to deduce the molecular events involved in the molting process. The results provided a basic understanding of the molecular mechanisms undergirding molting, including hormone regulation, skelemin, immune response etc. Our results show that molting cycles of L. vannamei are encoded by a large number of gene families subject to strict patterns of temporal and spatial regulation. We propose a cascade of molecular events applicable to the molting cycle of L. vannamei (Fig 8). These comprehensive analyses provide molecular evidence that will improve our understanding of morphological variation in molting and serve as a potential blueprint for future research on molting in crustaceans and other molting animals.  Table. Expression profiles of all 5,117 DEGs from the L. vannamei molting cycle. Each row is a DEG identified from the molting transcriptome. Expression levels in all eight stages are provided for each unigene. (XLSX) S4 Table. Details of 259 putative pathways identified for the eight stages of molting. The name of the pathway, its hierarchy, and the number of members plus their gene IDs are provided. (XLSX) S5 Table. Thirty-six immunity families, including 236 immune genes, categorized between molting stages. Gene ID, number of members and expression for the eight molting stages are given for each unigene. (XLSX)