Transcriptomic variation of eyestalk reveals the genes and biological processes associated with molting in Portunus trituberculatus.

BACKGROUND: Molting is an essential biological process throughout the life history of crustaceans, which is regulated by many neuropeptide hormones expressed in the eyestalk. To better understand the molting mechanism in Portunus trituberculatus, we used digital gene expression (DGE) to analyze single eyestalk samples during the molting cycle by high-throughput sequencing.
RESULTS: We obtained 14,387,942, 12,631,508 and 13,060,062 clean sequence reads from inter-molt (InM), pre-molt (PrM) and post-molt (PoM) cDNA libraries, respectively. A total of 1,394 molt-related differentially expressed genes (DEGs) were identified. GO and KEGG enrichment analysis identified some important processes and pathways with key roles in molting regulation, such as chitin metabolism, peptidase inhibitor activity, and the ribosome. We first observed a pattern associated with the neuromodulator-related pathways during the molting cycle, which were up-regulated in PrM and down-regulated in PoM. Four categories of important molting-related transcripts were clustered and most of them had similar expression patterns, which suggests that there is a connection between these genes throughout the molt cycle.
CONCLUSION: Our work is the first molt-related investigation of P. trituberculatus focusing on the eyestalk at the whole transcriptome level. Together, our results, including DEGs, identification of molting-related biological processes and pathways, and observed expression patterns of important genes, provide a novel insight into the function of the eyestalk in molting regulation.

Introduction

Molting is an essential biological process occurring multiple times throughout the life history of crustaceans, and is essential for development, growth, and reproduction [1, 2]. Crustaceans experience rhythmic molting cycles that contain three main stages, inter-molt (InM), pre-molt (PrM), and post-molt (PoM), that are based on morphological and microstructure observations. The InM stage involves energy accumulation and muscle regeneration. During the PrM stage, the old exoskeleton is reabsorbed and a new exoskeleton is formed. In the PoM stage, water uptake leads to expansion, which is essential for growth.

A possible signaling pathway in molting regulation has been proposed, which includes triggering and summation phases, and involves the decapod crustacean molting gland [3]. Many regulatory neuropeptides and neuropeptide receptors and other crucial genes pertaining to molting have been identified, including molt-inhibiting hormone (MIH), crustacean hyperglycemic hormone (CHH), retinoid X receptor and ecdysteroid receptor (EcR) [4-8], cuticle-related enzymes (chitinases, chitin deacetylase, and chitin synthase) [9-11], and structural proteins of the cuticle [12, 13]. Although many studies focusing on crustacean molting have been carried out, our knowledge of the molecular mechanisms involved is still far from complete.

The eyestalk of crustaceans, in which the X-Organ/Sinus Gland complex (XO-SG) is located, plays important roles in crustacean molting [14]. The XO-SG complex releases a number of peptide hormones into the hemolymph, and these hormones are involved in a range of physiological activities such as molting, growth and postembryonic development in crustaceans [15]. The MIH triggers molting upstream of a signaling pathway that links MIH to the regulation of ecdysteroidogenesis in the decapod crustacean molting gland [16], which is secreted by the XO-SG complex [3, 17]. The multifunctional neuroendocrine hormones and CHH, which are also synthesized and stored in the XO-SG complex, were found to inhibit crustacean molting [18]. The crustacean eyestalk is covered by the extracellular cuticle, an important downstream target effector organ, which is thus an ideal organ to study transcriptomic changes during molting cycles.

In recent years, in order to clarify the molecular mechanisms of molting regulation, some molting-related genes in crustaceans have been identified via microarray technology [12, 13]. However, due to the limitations of the technology, hybridization probes can only be designed for limited cDNA or expressed sequence tag data, leading to a huge number of unknown genes or transcripts that cannot be detected [19]. RNA-seq via high-throughput sequencing is a powerful technology for identifying genes and has several distinct advantages, such as high throughput, higher sensitivity and the ability to detect new or unknown genes [13]. Due to the importance of the eyestalk in crustaceans, significant progress has been made in understanding the profile of transcript expression via RNA-seq over the last 2 years; however, these studies have mainly focused on development and neuropeptide mining rather than clarifying the mechanism of molting regulation [20-22].

Portunus trituberculatus (Crustacea: Decapoda: Brachyura), commonly known as the swimming crab, is widely distributed in the coastal waters of China, Japan, Korea and other East Asian countries [19]. The crab has become one of the most important economic species for its high nutritional value and fast growth in marine aquaculture [23]. Molt stages of the crab are divided into four basic periods depending on morphological feature observation [24]: post-molt (stage A and stage B), inter-molt (stage C), pre-molt (substage D0, substage D1, substage D2, substage D3, substage D4) and molt (stage E). Young juvenile crabs experience 8–10 moltings before sexual maturity [25], and show saltatory growth after each molting cycle. Comprehensive understanding of the molecular mechanisms of molting regulation is essential for controlling its growth and development of P. trituberculatus, which would benefit the aquaculture industry in China.

To investigate and gain a better understanding of the molting mechanism in P. trituberculatus, we created a reference transcriptome and digital gene expression (DGE) profile for single eyestalk samples in three major molting stages (InM, PrM, and PoM) using Illumina sequencing. To our knowledge, this is the first transcriptome-wide gene expression profiling of eyestalk in P. trituberculatus. This study will be a foundational resource for further studies in molting regulation mechanisms.

Methods

Animals

Ninety healthy swimming crabs (P. trituberculatus) between 80 and 100 days of age were obtained from Haifeng Company (Changyi, Weifang, China). The crabs used in the present study were artificial breeding animals. This species is the most common edible crustacean in China and is not an endangered or protected species. All the samples were acclimated in the laboratory (33 ppt, 18°C) for 1 week before beginning the experiment. Three molt stages were identified based on morphological features [24] and the expression level of MIH mRNA [26], which included intermolt (InM, stage C, all parts of the body were hard and the new carapace had not yet started secretion), premolt (PrM, substage D3/D4, the old and new carapaces were separated completely) and postmolt (PoM, stage A, parts of the body were flaccid). Nine male crabs from each molt stage (InM, PrM, and PoM) were selected. Then the crabs were placed in an ice bath until anesthetization (5–10 min), and eyestalk were dissected, cleaned of pigments and exoskeleton, and then stored in RNAlater (−20°C) until RNA extraction.

RNA extraction, library preparation, and sequencing

Total RNA was extracted individually using Trizol (Invitrogen, Carlsbad, CA, USA) and then equivalent RNAs were pooled into each group. For library preparation, 3 μg RNA per group was used as input material. Sequencing libraries were generated using a NEBNext Ultra RNA Library Prep Kit for Illumina (NEB, USA) following the manufacturer’s recommendations and index codes were added to attribute sequences to each sample. Library quality was assessed on the Agilent Bioanalyzer 2100 system. The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumina) according to the manufacturer’s instructions. Then, the library preparations were sequenced on an Illumina Hiseq 2000 platform and 100 bp single-end reads were generated.

Assembly and annotation

The raw reads from our two previous works were used to assemble the reference transcriptome [19, 27]. The clean reads were assembled by Trinity as described previously [28], followed by TIGR Gene Indices clustering tools (TGICL) [29]. The longest assembled sequences were referred to as contigs. The reads were then mapped back to contigs with paired-end reads to detect contigs from the same transcript and the distances between these contigs. Finally, sequences were obtained that lacked Ns and could not be extended on either end [30]. Such sequences were defined as unigenes. The unigenes were annotated via public databases (Nr; Nt; Swiss-Prot; COG; and Kyoto Encyclopedia of Genes and Genomes, KEGG) by BlastX via an E-value cut-off of 1.0 x 10−5.

DEG identification

The read counts were obtained by RSEM after being mapped back to the reference transcriptome, which was normalized via RPKM [30]. Prior to differential gene expression analysis, for each sequenced library, the read counts were adjusted by the edgeR program package through one scaling normalized factor [31]. Differential expression analysis of two conditions was performed using the DEGSeq R package (1.12.0) [32]. The P values were adjusted using the Benjamini & Hochberg method. Corrected p-value of 0.005 and log2 (fold change) of 1 were set as the threshold for significantly differential expression. Selected differentially expressed genes (DEGs) were clustered by the STEM Clustering Method [33].

GO and KEGG pathway enrichment analysis

Gene Ontology (GO) enrichment analysis of DEGs was implemented by the GOseq R package, in which gene length bias was corrected. GO terms with a corrected P value of less than 0.05 were considered significantly enriched for DEGs. KEGG is a database resource for understanding high-level functions and utilities of the biological system, such as the cell, the organism, and the ecosystem, from molecular-level information, especially large-scale molecular datasets generated by genome sequencing and other high-throughput experimental technologies (http://www.genome.jp/kegg/). We used KOBAS to test the statistical enrichment of differential expression genes in KEGG pathways [34].

Real-time PCR

To test the reliability of RNA-seq, 18 unigenes were selected to confirm the accuracy of different expression patterns. Their specific primers were designed with Primer Premier 5 software (Premier Biosoft International) (S1 Table). First strand cDNA was synthesized with prepared RNA of Illumina sequencing using a PrimeScript RT reagent kit (Takara, Dalian, China). Quantitative real-time PCR was performed using a 7500 Fast Real-Time PCR System and QuantiFast SYBR Green PCR Kit (Qiagen, USA). The PCR was carried out in a total volume of 10 μl and performed with the following thermal profile: 95°C for 5 min, 40 cycles of 95°C for 10 s and 60°C–65°C for 30 s. Fluorescence levels were measured after the 60°C step. At the end of the PCR cycles, a melt curve was generated to analyze product specificity.

Results

Reference transcriptome assembly and annotation

In order to achieve a comprehensive P. trituberculatus reference transcriptome, previous RefSeq data obtained from a variety of tissues (including the eyestalk, gill, heart, hepatopancreas, and muscle) by Illumina deep sequencing was used to create the Trinity-assembled transcriptome. The overall Illumina sequencing data were deposited in the Short Read Archive database of the National Center for Biotechnology Information (NCBI) (SRR1013694, SRR1013695, SRR1013696, SRR1168416, and SRR1168417), which contained 155,713,331 read pairs and 31.02G clean bases (Table 1). Finally, we identified 186,081 transcripts with an average length of 1,028 bp (Fig 1), which contained 130,560 unigene sequences (submitted in TSA database of NCBI: SUB2309547).

Table 1

Summary of Illumina transcriptome sequencing, assembly and annotation for Portunus trituberculatus.

Raw results (after trimming)		Assembly results		Annotation results
Clean bases (G)	31.02	Transcripts	186,081	Nr annotations	17,639
Read pairs	155,713,331	Average length of transcripts (bp)	1,028	KOG hits	16,010
Read length (bp)	100	Smallest transcripts (bp)	201	GO mapped	33,723
		Largest transcripts (bp)	36,343	KEGG hits	8,087

Fig 1

Sequence length distribution of transcripts assembled from Illumina reads.

Overall, 33,723 (25.83%) unigene sequences were categorized into 58 subcategories according to GO [35]. Six subcategories, ‘cellular process’, ‘binding’, ‘cell’, ‘cell part’, ‘metabolic process’ and ‘catalytic activity’ were dominant clusters in the GO classification (Table 1; S1 Fig). A total of 16,010 unigenes could be grouped into 26 functional classes via the EuKaryotic Ortholog Groups program (KOG) [36]. Among which, ‘general function prediction’ (2,926), ‘signal transduction’ (2,488), and ‘post-translational modification, protein turnover, chaperone’ (1,720) were the dominant groups and represented 18.28%, 15.54%, and 10.74%, respectively (Table 1; S2 Fig); 8,087 unigenes were classified to 38 KEGG pathways using the KEGG database [37]. ‘infectious diseases’ (1,619 members) and ‘translation’ (1,075 members) were the two largest pathways (Table 1; S3 Fig).

DEG identification and analysis during different molting cycles

To identify DEGs for eyestalks in different molting cycles, we sequenced three eyestalk libraries for each molting stage (InM, PrM, and PoM). Finally, the three libraries generated 14,387,942, 12,631,508 and 13,060,062 clean reads, which were mapped to the reference transcriptome and represented 74.20%, 74.10%, and 77.49%, respectively (S2 Table). The overall Illumina sequencing data were deposited in the Short Read Archive database of the NCBI (SRR5166969, SRR5166967 and SRR5166965).

A total of 1,394 DEGs were identified by DESeq [38]. The number of genes that were up- or down-regulated at the different molt stages is shown in Fig 2. There were 445 DEGs during the InM and PrM stages, among which, 284 were down-regulated and 161 were up-regulated in PrM stages, and 1,181 DEGs (452 down-regulated and 666 up-regulated) were identified between the PrM and PoM stages. A total of 728 DEGs (476 down-regulated and 252 up-regulated) were identified between InM and PoM.

Fig 2

Venn diagrams of the differentially expressed genes among three molting stages.

Eighteen unigenes were subjected to qPCR assays in order to verify the results of differential gene expression analysis (Table 2). Within the 54 comparisons, 40 (74.1%) matched well with the DEG data. Although qPCR data of 14 unigenes did not match the DEG data, they showed similar trends with the Illumina sequencing in up- or down-regulation. The results validated that the majority of the DEGs identified from the DEG analysis in different molting phases were authentic.

Table 2

Validation of DEGs data for the selected eighteen transcripts using qPCR for relative gene expression.

Gene ID	Gene Name	Illumina sequence			qPCR
		InM vs. PrM	PoM vs. PrM	PoM vs. InM	InM vs. PrM	PoM vs. PrM	PoM vs. InM
comp101914_c0	hemocyanin beta subunit 1	1.27	2.28	0.78	3.67*	2.14	0.64
comp99362_c0	cAMP-responsive element-binding protein-like 2-like	1.01	0.53	0.52	1.17	0.35	0.39
comp103008_c0	neurogenic locus Notch protein-like	0.61	0.47	0.76	0.25*	0.34	0.72
comp98266_c0	transferrin	0.64	0.52	0.80	0.54	0.58	0.59
comp96475_c1	Espin	0.82	13.37*	16.30*	0.43	9.94*	23.09*
comp102262_c0	calpain-c	0.70	8.73*	12.32*	0.86	9.2*	8.03*
comp99514_c0	transient-receptor-potential-like protein	0.48	0.01*	0.01*	0.41	0.01*	0.01*
comp102473_c0	Sodium/calcium exchanger 1	0.70	1.33	1.89	0.94	2.03	1.40
comp57954_c0	cuticle protein VER3-like	0.69	4.77*	6.84*	0.66	10.03*	5.90*
comp85008_c0	Chitin binding Peritrophin-A domain	0.84	6.40*	7.03*	0.50	11.09*	3.13*
comp90850_c0	cuticle protein CB2-like	0.31*	5981.73*	18839.83*	0.19*	8172.91*	13834.76*
comp88516_c0	cuticle protein CB1	0.36	6326.38*	17268.64*	0.49	1656.17*	12670.09*
comp61465_c0	cuticle protein	0.63	8665.93*	13646.97*	0.65	1802.63*	13425.21*
comp36475_c0	cuticle protein BD1	2.38	11898.15*	4996.53*	5.06*	2644.29*	4940.72*
comp90588_c0	HR4 nuclear receptor	0.65	3.20*	4.93*	0.37	1.48	4.02*
comp83668_c0	Probable nuclear hormone receptor HR38	1.59	3.88*	2.43	1.23	1.09	1.88
comp83095_c1	Nuclear hormone receptor FTZ-F1	1.35	6.70*	4.94*	1.65	1.94	2.99
comp93173_c0	Ecdysone-induced protein 74EF isoform B	0.73	0.48	0.65	0.39	0.13*	0.35

* represents a significant difference between the InM, PrM and PoM.

GO and KEGG enrichment analysis

DEGs were further analyzed according to GO functional enrichment analysis (Fig 3). The 188 DEGs of PrM vs. InM were enriched in 10 processes, among which, most of the DEGs (77.66%) were up-regulated in PrM. In the 10 processes, chitin metabolic process (GO:0006030), ribosome (GO:0005840) and peptidase inhibitor activity (GO:0030414) have the maximum enrichment level in biological process (BP), cellular component (CC), and molecular function (MF), respectively. The 393 DEGs of PoM vs. PrM were enriched in nine processes, which included 157 down-regulated and 236 up-regulated DEGs in PoM, respectively. From PrM to PoM, chitin metabolic process (GO:0006030, BP), extracellular region (GO:0005576, CC) and structural constituent of cuticle (GO:0042302, MF) were processes with the highest enrichment level. The 268 DEGs of InM vs. PoM were enriched in seven processes, among which, most of the DEGs (79.10%) were down-regulated in InM. Chitin metabolic process (GO:0006030, BP), extracellular region (GO:0005576, CC) and structural constituent of cuticle (GO:0042302, MF) were processes with the highest enrichment level between PrM and PoM.

Fig 3

GO enrichment analysis of the differentially expressed genes.

DEGs were also analyzed according to KEGG enrichment analysis in a complete molting cycle. Six groups of DEGs were enriched in 34 KEGG pathways (Fig 4), of which DEGs up-regulated in PrM vs. InM were enriched to the largest number of KEGG pathways (18). During the PoM and PrM stages, the up- and down-regulated DEGs were enriched in six and 15 KEGG pathways, respectively. The pathway of 'amino sugar and nucleotide sugar metabolism' (ko00520) was enriched at the highest level (p-value = 2.43×10−9). Eighteen KEGG pathways were enriched during PoM and InM stages, among which, 13 KEGG pathways were up-regulated and five KEGG pathways were down-regulated.

Fig 4

KEGG enrichment analysis of the differentially expressed genes.

Group 1, DEGs up-regulated in PrM vs. InM. Group 2, DEGs down-regulated in PrM vs. InM. Group 3, DEGs up-regulated in PoM vs. PrM. Group 4, DEGs down-regulated in PoM vs. PrM. Group 5, DEGs up-regulated in InM vs. PoM. Group 6, DEGs down-regulated in InM vs. PoM.

Expression pattern analysis of molting-related gene

Forty-seven molting-related genes were screened and clustered by the STEM Clustering Method [33], which included 6 neuropeptide transcripts, 4 molting-related receptor transcripts, 31 cuticle transcripts and 6 chitinase transcripts. (Fig 5 and Table 3). They were grouped into four main clusters (labeled 1–4) according to their expression patterns via K-means clustering, among which, Profile 1 contained the highest number of transcripts (36, 76.6%), followed by Profile 2 (9), Profile 3 (1) and Profile 4 (1). Five neuropeptide transcripts (83.3%) were clustered in Profile 2, which were up-regulated in PrM, and then down-regulated in PoM. Most of the cuticle transcripts (29, 93.5%) and chitinase transcripts (5, 83.3%) were clustered in Profile 1, which were not significantly changed in PrM to InM and up-regulated in PoM. Four molting-related receptor transcripts showed different expression patterns, which were clustered in Profile 1 (2) and 2 (2).

Fig 5

Cluster analysis of the 47 molting related genes.

Forty-seven genes were grouped into four main clusters (labelled Profile 1–4) according to their expression patterns via K-means clustering. The three dots in order represented InM, PrM and PoM respectively. Lines in different colours represented different genes.

Table 3

Information on the 47 clustered genes (DEGs).

Gene_id	Gene categories	Profile
comp103393_c0	cuticle	Profile 1
comp104448_c0	cuticle	Profile 1
comp57738_c0	cuticle	Profile 1
comp94993_c0	cuticle	Profile 1
comp79334_c1	cuticle	Profile 1
comp89919_c0	cuticle	Profile 1
comp62755_c0	cuticle	Profile 1
comp99348_c0	cuticle	Profile 1
comp95871_c1	chitinase	Profile 1
comp36064_c0	cuticle	Profile 1
comp95871_c0	chitinase	Profile 1
comp89113_c1	cuticle	Profile 1
comp108916_c0	chitinase	Profile 1
comp83095_c1	molting related receptor	Profile 1
comp53926_c0	cuticle	Profile 1
comp36475_c0	cuticle	Profile 1
comp35636_c0	cuticle	Profile 1
comp101437_c0	cuticle	Profile 1
comp138619_c0	cuticle	Profile 1
comp302044_c0	cuticle	Profile 1
comp81349_c1	cuticle	Profile 1
comp81691_c0	cuticle	Profile 1
comp70997_c0	cuticle	Profile 1
comp91479_c0	cuticle	Profile 1
comp91757_c0	cuticle	Profile 1
comp86312_c0	cuticle	Profile 1
comp61465_c0	cuticle	Profile 1
comp58768_c0	cuticle	Profile 1
comp88516_c0	cuticle	Profile 1
comp78778_c0	chitinase	Profile 1
comp96722_c0	cuticle	Profile 1
comp99034_c0	chitinase	Profile 1
comp63340_c0	cuticle	Profile 1
comp90594_c0	cuticle	Profile 1
comp90588_c0	molting related receptor	Profile 1
comp104537_c0	cuticle	Profile 1
comp100803_c0	neuropeptide	Profile 2
comp103698_c3	neuropeptide	Profile 2
comp92243_c0	neuropeptide	Profile 2
comp96149_c0	neuropeptide	Profile 2
comp97030_c0	neuropeptide	Profile 2
comp86620_c0	cuticle	Profile 2
comp103455_c0	chitinase	Profile 2
comp86338_c0	molting related receptor	Profile 2
comp97869_c0	molting related receptor	Profile 2
comp90990_c0	neuropeptide	Profile 3
comp37103_c0	cuticle	Profile 4

Discussion

Molting is a vitally important metamorphosis process for P. trituberculatus and affects its physical growth, tissue regeneration, and gonad development. Apart from studies concentrating on MIH, molt hormone, and methyl farnesoate, there has been little molecular research on the molting process [4]. Our study is the first to use high-throughput filtration to investigate molting-related genes in the eyestalk of P. trituberculatus at the whole transcriptome level. Compared with other technologies, RNA-seq technology has several advantages, such as high throughput, higher sensitivity and the ability to detect new or unknown genes. This approach will assist in the discovery and annotation of novel genes that play key roles in the molting process.

The reference transcriptome used in this study was based on five groups of transcriptome data from our two previous works [19, 27], in which equal quantities of RNA were extracted and mixed from five tissues, including the eyestalk, and sequenced by Illumina. The transcript library obtained covered all of the P. trituberculatus ESTs from NCBI [27]. The overall Illumina read pairs and clean bases for all samples used in this study are 155,713,331 and 31.02G. We finally identified 186,081 transcripts after assembly analysis, the number of transcripts was significantly higher compared to the previous studies (120,137 and 141,339, respectively). Therefore, the transcriptome assembled in this study could be used as a reference transcriptome with a high coverage rate for subsequent analysis.

Identification of DEGs among the different molting phases is important to find the underlying candidate genes of molting regulation in P. trituberculatus. We detected a total of 1,394 DEGs when comparing the three important molting cycles (qvalue<0.005 & |log2 (fold change)|>1). Some genes with important roles in the regulation of molting were found in the DEGs, which included an ecdysteroid receptor (comp97869_c0) [39], chitinase (comp99034_c0 and comp108916_c0) [9, 40], hemocyanin (comp101914_c0) and nuclear hormone receptor (comp86338_c0 and comp86338_c0) [41]. Many of the DEGs have not yet been reported to have functions in the molting cycle, in spite of their having important roles in other physiological processes. For example, heat shock protein 70 (comp84939_c0), sodium/hydrogen exchanger (comp101931_c0) and methyltransferase (comp83893_c0) play an important role in immunity, osmoregulation, and epigenetic modification, respectively. This suggests that these genes might have multiple functions and may also play important roles in the regulation of molting. It also suggests that the process of molting is closely related to these physiological processes. It should be noted that most of the DEGs (746/1394) could not be aligned to any known protein by BlastX (E values less than 10−5). This indicates that the regulation of molting was far more complicated than we originally thought, and that there is still a long way to go in clarifying the molting mechanism in crustaceans.

GO and KEGG enrichment analysis were able to improve our understanding of the underlying biological processes related to various biological traits. A total of 15 GO processes were enriched, of which five processes (chitin metabolic process, peptidase inhibitor activity, chitin binding, carbohydrate derivative binding and structural constituent of cuticle) were found in every comparison of the three groups, which suggests these processes play an important role in molting. Thirty-four pathways were enriched via KEGG enrichment analysis. From PoM to PrM, the enriched pathways were mainly up-regulated (85.7% and 72.2% in InM vs. PoM and PrM vs. InM, respectively), which included ‘ribosome’, ‘amino sugar and nucleotide sugar metabolism’, ‘retrograde endocannabinoid signaling’, ‘circadian entrainment’ and ‘phosphatidylinositol signaling system’. The pathway of ‘ribosome’, with the highest enrichment level and the highest number of DEGs (p-value = 2.84461E-10, 24 DEGs), was enriched in PrM vs. InM. The ribosome is a complex molecular machine, found within all living cells, that serves as the site of biological protein synthesis (translation). Our results indicate that a large number of protein synthesis processes were triggered to prepare for molting during PrM. From PrM to PoM, 71.4% KEGG pathways were down-regulated, which suggested that the body was in a calm state after molting. However, the pathway of ‘amino sugar and nucleotide sugar metabolism’ with the highest enriched level (p-value = 4.14001E-08) was up-regulated in PoM vs. PrM. An amino sugar is a sugar molecule in which a hydroxyl group has been replaced with an amine group. More than 60 amino sugars are known, with one of the most abundant being N-acetyl-d-glucosamine, which is the main component of chitin. This indicated that there could be more active cuticle metabolism in the body after the early phase of PoM.

Previous studies on Eriocheir sinensis found a DEGs pattern associated with energy metabolism during a molting cycle. To be specific, DEGs enriched in PoM were linked to energy consumption, whereas genes enriched in InM were related to carbohydrates, lipids metabolic and biosynthetic processes [42]. However, we did not find a similar pattern with energy metabolism in our study. We believe this difference is mainly caused by the different tissues used in the two studies. Unlike the hepatopancreas, which has a main role in energy metabolism and storage, the eyestalk is a major site for the production of neurohormones and controls a variety of physiological functions such as molting and reproduction [20]. Therefore, it was not surprising that there were little energy metabolism processes or pathways enriched in DEGs in this study. But it is worth noting that a number of neuromodulator-related pathways were enriched including the dopaminergic, serotonergic, cholinergic and glutamatergic pathways, which were not found in previous research [12, 13, 42]. Hormones from the XO-SG complex appear to be controlled by the neurotransmitters/neuromodulators dopamine, serotonin, and enkephalin [43-47]. The release of some SG hormones is controlled presynaptically by dopamine, metenkephalin, and leuenkephalin, which occurs via synaptic contacts in the medulla terminalis near neuron cell bodies [48]. These neuromodulator-related pathways were up-regulated in PrM and down-regulated in PoM, which indicates molting-related neurohormone release is triggered by neuromodulators in the early phase of PrM, then aroused a series of subsequent molting processes. However, the relationship between neurotransmitters and molting-related hormones is far form clear, and further research is necessary to clarify the underlying mechanisms regulating hormone release.

The process of molting is regulated by an elaborate interplay of neuropeptide hormones, which include positive and negative regulation [13]. The signals that arise by neuropeptide hormones trigger a series of physiological activities related to molting, including losing the extracellular cuticle, escaping from the confines of the cuticle relatively rapidly, taking up water, expanding the new flexible exoskeleton and then quickly hardening it for defense and locomotion. In particular, meticulous control of regulatory proteins and genes is required to form and harden a new exoskeleton as well as to dissolve, reabsorb and shed the old one [2, 49]. Forty-seven DEGs were screened and clustered, which included six neuropeptide transcripts, four nuclear receptor (NR) transcripts, 31 cuticle transcripts and six chitinase transcripts. Most genes had similar expression patterns (Profile 1 and 2) via clustering analysis. This suggests that there is a connection between these genes throughout the molt cycle.

Crustacean neuropeptides and NRs play important roles in regulating many physiological activities, including molting, growth, feeding and reproduction [50-52]. In arthropods, the molting process is regulated by neuropeptide hormones acting via NR proteins [41]. Our results show that the highest number of these two kinds of genes were clustered in profile 2, including five neuropeptides (pigment-dispersing hormone [PDH], neuroparsins [NPs] 1 and 2, neuropeptide and neuroendocrine convertase) and 2 NRs (EcR and E75). Interestingly, the five neuropeptide transcripts were significantly up-regulated in PrM, and their functions in molting regulation have rarely been reported. For example, the PDH, which is released primarily from the eyestalk, is important in regulating pigment dispersion in crustaceans [53]. NPs were initially identified in insects with diverse roles including anti-juvenile, antidiuretic and hyperlipidemic effects [54, 55]. Neuropeptide F (NPF) is a potent stimulator regulating feeding behavior, circadian rhythm and reproduction [56-59]. The five neuropeptides were up-regulated in PrM vs. InM and subsequently down-regulated in PoM vs. PrM indicating that they play important roles in regulating molting in PrM.

EcR belongs to the NR superfamily and acts as an important transcription factor regulating downstream genes involved in the molting process [60]. Consistenting with previous studies in crustacean molting [61, 62], the expression of EcR in this research was significantly up-regulated in PrM and then down-regulated in PoM. E75, which also belongs to the NR superfamily, has been identified in numerous insects such as D. melanogaster, Aedes aegypti and Galleria mellonella [63-65]. Although the complex role of E75 is not fully understood, it was proved to be an responsive gene in ecdysteroids signaling pathway in crustacean [66, 67]. In this study, the expression pattern of E75 was similar to that of EcR, and both of them were clustered in profile 2. Our results suggest that EcR and E75 play essential roles in molting of P. trituberculatus, and point to a possible synergy between them, which deserves further verification.

Chitin, a 1, 4-β-linked polymer of N-acetyl-β-glucosamine, is found in crustacean shells [68]. During the pre-molt period, the epidermis secretes chitinases that degrade the inner layers of the old exoskeleton while synthesizing a new exoskeleton [69]. In this study, one of the six chitinase transcripts was clustered in Profile 1, which was up-regulated in PrM vs. InM and down-regulated in PoM vs. PrM, suggesting a specific role in degrading the old exoskeleton in the PrM stage. However, the chitinases have diverse functions in a wide range of organisms, including morphogenesis [70], host defense against chitin-bearing pathogens, cell division and sporulation [71], and regulators of development [72]. Interestingly, five of the six chitinase transcripts were clustered in Profile 1, which were significantly up-regulated in PoM, indicating they may participate in other important physiological processes after the molting stage, such as growth and development [73].

The crustacean cuticle provides initial reinforcement by cross-linking cuticular proteins attached to the cuticle chitin-fiber matrix and the protein matrices generally have very complex compositions. Thirty-one cuticle proteins were differentially expressed across the molt cycle and mainly clustered in Profile 1. The majority of them have a common variation tendency: they were not significantly changed in PrM vs. InM, but significantly up-regulated in the PoM, which shows that the majority of cuticles in P. trituberculatus were highly active at the PoM stage. The endocuticle and membranous layer form during the PoM, accompanied by calcification and sclerotization [74]. Our data suggest that the process to generate a new cuticle mainly occurs in the PoM, which is similar to other reports of crustacean cuticle protein expression during the molting cycle [69, 75].

Conclusions

Our work is the first molt-related investigation of P. trituberculatus that is focused on the eyestalk at the whole transcriptome level. A total of 1,394 molt-related DEGs were identified. Many important processes and pathways that play a key role in molting regulation were identified including chitin metabolism, peptidase inhibitor activity, ribosome and circadian entrainment. In particular, we were the first to observe a pattern associated with neuromodulator-related pathways and four kinds of important molting-related gene families during a molting cycle. Together, our results, including screened DEGs, identified molting-related biology processes and pathways, and the observed expression pattern of important genes, provide novel insights into the functions of the eyestalk in molting regulation.

GO classifications for the assembled unigenes.

(TIF)

Click here for additional data file.

KOG classifications for the assembled unigenes.

(TIF)

Click here for additional data file.

KEGG classifications for the assembled unigenes.

(TIF)

Click here for additional data file.

List of Differently Expressed Genes (DEGs).

(XLSX)

Click here for additional data file.

Summary of Digital Gene Expression (DGE).

(XLSX)

Click here for additional data file.