Comparative transcriptomic analysis provides insights into the molecular basis of the metamorphosis and nutrition metabolism change from zoeae to megalopae in Eriocheir sinensis

Comparative Biochemistry and Physiology, Part D 13 (2015) 1–9

Contents lists available at ScienceDirect

Comparative Biochemistry and Physiology, Part D journal homepage: www.elsevier.com/locate/cbpd

Comparative transcriptomic analysis provides insights into the molecular basis of the metamorphosis and nutrition metabolism change from zoeae to megalopae in Eriocheir sinensis Yingdong Li a,1, Min Hui a,1, Zhaoxia Cui a,b,⁎, Yuan Liu a, Chengwen Song a,c, Guohui Shi a,c a b c

EMBL, Institute of Oceanology, Chinese Academy of Sciences, Qingdao, 266071, China National and Local Joint Engineering Laboratory Of Ecological Mariculture, Qingdao, 266071, China Graduate University of the Chinese Academy of Sciences, Beijing, 100049, China

a r t i c l e

i n f o

Article history: Received 8 August 2014 Received in revised form 27 October 2014 Accepted 27 October 2014 Available online 4 November 2014 Keywords: Larval development Metamorphosis Mitten crab Nutrition metabolism Transcriptomics

a b s t r a c t Within the larval period of Eriocheir sinensis, there is pronounced morphological changes upon the molt from the fifth zoeae (Z5) to megalopae (M), and low survival rate exists during this transition, which is typical in crab species. RNA sequencing was applied to Z5 and M of E. sinensis, resulting in the discovery of 19,186 unigenes and 652 differentially expressed genes (DEGs, 3.40% of the unigenes). The important metabolic pathways that might play roles in the larval development of E. sinensis from Z5 to M were detected to be ‘Xenobiotics Biodegradation and Metabolism (8.16%)’, ‘Metabolism of Cofactors and Vitamins (6.70%)’, ‘Lipid Metabolism (6.36%)’, and ‘Amino Acid Metabolism (6.28%)’. Further, 19 DEGs possibly contributing to the morphological and sensory capability changes of the larvae were identified, like multiple copies of cuticle protein genes, retinaldehyde-binding protein 1 (RLBP1), envelope protein (Envelope) and hormone-related gene ecdysteroid-regulated 16 kDa protein (ESR16). Moreover, 62 DEGs were identified to be related to carbohydrate, lipid and protein digestion and metabolism, such as glucose dehydrogenases (GDHs), lipases (LIPs) and serine proteases (SPs). Among these DEGs, more genes related to the substance metabolism were found up-regulated in Z5 than M, suggesting that more energy might be essential to be released for Z5 to complete the transition into M. Characterization of the crucial DEGs by real-time quantitative PCR re-conformed their expression pattern. This study provides the first genomewide transcriptomic analysis of E. sinensis Z5 and M for studying the molecular basis of the larvae metamorphosis and nutrition metabolism. © 2014 Elsevier Inc. All rights reserved.

1. Introduction Studies of the early ontogenetic stages in benthic aquatic invertebrates are especially helpful for understanding the type of population control acting upon a species (Roughgarden et al., 1988; Osman and Whitlatch, 2004). This is due to the fact that the larvae and early juvenile stages are the most important for dispersal and they have the highest level of mortality at the population level (Palmer et al., 1996; Gosselin and Qian, 1997; Hunt and Scheibling, 1997). Additionally, in the case of species with commercial value, these studies allow us to understand the fluctuations in stock and recruitment (Wahle, 2003), which can supply useful information in the sustainable management of fisheries. Chinese mitten crab, Eriocheir sinensis (Crustacea, Decapoda, Grapsidae), is an important economic decapod crustacean species in

⁎ Corresponding author at: EMBL, Institute of Oceanology, Chinese Academy of Sciences, 7 Nanhai Road, Qingdao 266071, China. Tel./fax: +86 53282898509. E-mail address: [email protected] (Z. Cui). 1 These authors contributed equally to the study.

http://dx.doi.org/10.1016/j.cbd.2014.10.002 1744-117X/© 2014 Elsevier Inc. All rights reserved.

China and has spread to Europe and America as an invasive species (Herborg et al., 2003). The species is easily propagated and can be transported over long distance. Therefore, it is expected to serve it as a model species for decapod crustaceans study. The larval development of E. sinensis normally consists of five zoeal stages and a megalopal stage (Montú et al., 1996) as many other crab species, which correspond to the zoeal and mysis larval stages in shrimps. Two significant metamorphic processes exist in the development cycle of the crab, including molting from Z5 to M and from M to the first juveniles. When M are compared with Z5, the remarkable morphological features of M are the two apparent compound eyes, the stronger pereiopods and claws as well as the occurrence of pleopods after metamorphosis. Therefore, the megalopae have stronger activity capability, which can both swim and be benthic in float grass. Presently, knowledge regarding the molecular basis underling the metamorphosis in decapod crustaceans is scarce and mainly focused on exploring the process from larvae to the first juveniles instead of from zoeae to megalopae. In those studies, the 20-hydroxy ecdysone and methyl farnesoate (MF), a crustacean juvenile hormone, have been reported to play roles in molting and metamorphic transition (Laufer et al., 1987; Laufer and Biggers, 2001).

2

Y. Li et al. / Comparative Biochemistry and Physiology, Part D 13 (2015) 1–9

Recently, more genes and specific pathways have been revealed to be involved in the molting cycle and periodic morphological changes of the crustaceans using next-generation RNA sequencing technology (RNA-Seq) (Kuballa et al., 2011; Ventura et al., 2013). In the development of E. sinensis, low survival rate at metamorphosis from the fifth zoeal stage to the megalopal stage was notably found and has become a major problem in E. sinensis larvae culture (Cheng et al., 2008; Sui et al., 2011). It is believed that the low survival rate is principally resulted from the cannibalism in the larvae. The Z5 are usually predated by the newly metamorphosed M and the stronger individuals of megalopae can kill the weaker megalopae as well, which is mostly caused by asynchronous metamorphosis. In practice, the cannibalism and molting synchrony can be partially reduced by adjusting feeding ration and optimizing the feeding regime (Sui et al., 2011). Therefore, knowledge of the digestive physiology in the larvae is important for studies on the nutritional requirements and feeding ecology of crustacean (Biesiot and Capuzzo, 1990). However, most previous studies of food digestion in decapod crustaceans have been restricted to individual digestive enzyme activity (Lovett and Felder, 1990; Lemos et al., 1999; Johnston, 2003). Recently, by using RNA-Seq, potential genes involved in the synthesis and secretion of substances in another two crustacean species were further discovered based on transcriptomic sequences (Gonçalves et al., 2014; Wang et al., 2014). Studies on transcriptomes of E. sinensis have been performed in different adult tissues (Jiang et al., 2009; Zhang et al., 2011; He et al., 2012, 2013; Li et al., 2013), and only one transcriptome for the first zoeae at larval stage was reported (Cui et al., 2013). In our study, comparative transcriptomic analysis between the fifth zoeae and megalopae of E. sinensis was performed. Patterns of differential gene expression at the two developmental stages were screened in order to reveal the molecular basis of the principal change in the larvae at the two stages. This study also supplied more useful genetic resources for development biology study of crabs.

purified with AMPureXP beads and amplified for the construction of cDNA libraries. After cDNA libraries were prepared, the libraries were sequenced using the Illumina HiSeq™ 2000. 2.4. Transcriptome assembly Clean reads were obtained from raw reads by filtering the adaptor sequences, low-quality regions (bQ20) and sequences shorter than 50 bp using Solexa QA (Cox et al., 2010). In order to get more complete reference sequences, all clean reads were de novo assembled together with Trinity (http://trinityrnaseq.sourceforge.net/) (Grabherr et al., 2011). Trinity combines three independent software modules, Inchworm, Chrysalis and Butterfly, which were applied sequentially to process large volumes of RNA-seq reads into contigs, de Bruijin graphs and full-length transcripts. 2.5. Gene annotation and bioinformatics analysis Gene annotation was performed by sequence comparison with public databases. The transcripts were annotated using BlastX searches (cutoff E-value of 1E-05) against the non-redundant (NR) database (http://www.ncbi.nlm.nih.gov/) and then clustered into unigenes according to the top hit results. Then the Blast2GO program was used to obtain Gene Ontology (GO) terms and functional classification of all unigenes (Conesa et al., 2005). EC (Enzyme Commission number) terms and biochemical pathway information was generated by using Kyoto Encyclopedia of Genes and Genomes (KEGG) (http://www. genome.jp/kegg/) (Kanehisa et al., 2004). Evolutionary genealogy of genes of non-supervised orthologous groups (eggNOG) (http:// eggnog.embl.de/) was also performed to predict and classify potential functions of the unigenes based on known orthologous gene products (Powell et al., 2012). 2.6. Differential gene expression pattern analysis

2. Materials and methods 2.1. Ethical procedures The sampling location is not privately owned or protected, and no specific permission is required. No endangered or protected species were involved in the study. The experiments were performed in strict accordance with the guidelines set by the Institutional Animal Care and Use Committee (IACUC) of the Chinese Academy of Sciences (No. 2011-2). This study was specifically approved by the Committee on the Ethics of Animal Experiments of the Institute of Oceanology at the Chinese Academy of Sciences. All efforts were made to minimize the suffering of the larvae. 2.2. Preparation of experimental larvae Experimental larvae were collected from a commercial crab farm (Panjin Guanghe Crab Industry Co. Ltd.) in Liaoning, China, in June 2013. The larvae were checked by microscope to ensure their developmental stages. The larvae at Z5 and M stages were taken separately and immediately frozen in the liquid nitrogen. They were kept at − 80 °C for further use. Two replicated samples were prepared for Z5 and M, respectively. Thirty individuals were included in each sample. 2.3. cDNA library construction and deep sequencing Total RNA was extracted from the larvae of the samples according to the instructions of the Trizol Kit (Invitrogen, USA). The RNA was purified and cut into 155 bp fragments using TruSeq RNA Sample Prep Kit (illumina) following the manufacturer's protocol. Double-stranded cDNA was synthesized and sequencing adaptors were ligated as Illumina manufacturer’s instructions. The ligated products were then

The reads of Z5 and M samples were mapped back to the de novo assembly results separately, and gene expression profiling analysis was based on the number of the mapped reads. RPKM (reads per kb of exon model per million mapped reads) were used to calculate the expressed value (Mortazavi et al., 2008). DESeq program (http:// www-huber.embl.de/users/anders/DESeq/) was used to detect the differentially expressed genes (DEGs) between the two samples (Anders and Huber, 2010). DEGs were screened out with the standard expression fold change (N2) and significant adjusted p-value (p b 0.05), which was adjusted for the false discovery rate (FDR) due to multiple testing procedures to control the type I error (Benjamini and Hochberg, 1995). The DEGs were then submitted to GO enrichment analysis, KEGG Orthology (KO) and KEGG pathway enrichment analysis for further interpretation. According to the information from NR, GO and KEGG, crucial DEGs related to morphological change, hormone, nervous system and nutrition metabolism were further manually checked by searching public database and literatures according to the annotation from NR. 2.7. Real-time quantitative PCR (qRT-PCR) verification After total RNA from independent samples of Z5 and M was extracted separately, the first-strand cDNA was synthesised using M-MLV reverse transcriptase (Promega) and oligodT with 2 mg of total RNA. The cDNA was diluted 100 times by DEPC-treated water for the next step. The SYBR Green RT PCR assay was carried out in an ABI PRISM 7300 Sequence Detection System (Applied Biosystems). Six pairs of gene-specific primers (Mucin, RLBP1, PLA2, AGMO, PRSS and SP22D; Table 1) were used to amplify the partial cDNA gene sequences, respectively. Three biological replicates for each sample and three technical replicates were performed. The relative expression level was calculated

Y. Li et al. / Comparative Biochemistry and Physiology, Part D 13 (2015) 1–9 Table 1 Primer information for the qRT-PCR. Gene name

Unigene ID

Primer name

Primer sequences

Mucin

comp44879_c0_seq1

RLBP1

comp39285_c0_seq1

PLA2

comp47494_c4_seq7

AGMO

comp31174_c0_seq1

PRSS

comp44886_c0_seq1

SP22D

comp35785_c1_seq1

Mucin-F Mucin-R RLBP1-F RLBP1-R PLA2-F PLA2-R AGMO-F AGMO-R PRSS-F PRSS-R SP22D-F SP22D-R

GCATCACTAACTAATCGCCAC TCACACAGCGAACTGGTAAGG CGACCCAACCATTGACCAGC AGCACTGGAAGGCTTTGGAC ACTCTCCTACAAGGGCTACG TTCAGGCACTCCACGAACTC TACGCCGTACCACCTCAAGAG GCAGAAGTCCACGCCTATCAC CCACAGCAGTTCACGGTTCG CCGACGGTACTCCTCTGAAG CTCGCACTTAGCAGGCAAATG TCACAGCTGCTATGAGCCAAG

using the 2−ΔΔCt method. The β-actin gene was used to normalise the gene expression (Shen et al., 2014). The results were subjected to oneway analysis of variance (one-way ANOVA) using SPSS 16.0, and the p-values less than 0.05 were considered statistically significant. 3. Results and discussion 3.1. Phenotypic characterisation of the fifth zoeae and megalopae The larvae at Z5 and M stages were identified and characterised under the microscope according to those described by Montú et al. (1996). In Z5, the number of setae on the exopod of the maxilliped 1 and 2 was 12, while the number of setae at the telson was 10. These were used to distinguish Z5 from the other zoeal stages. When M were compared with Z5, the most strike features in M were the large stalked eyes and the presence of stronger pereiopods. Therefore, the M are relatively strong swimmers and walkers, and can prey on food actively with more acute sensory capability.

3

201 bp to 18,964 bp (Supplementary file 1) with an average length of 1,697 bp and an N50 of 2,629 bp (Table 2). Length distribution of the assembled contigs, transcripts and unigenes were displayed in Supplementary file 1. All raw reads sequence data from E. sinensis were deposited in the Sequence Read Archive (SRA) database (http://www. ncbi.nlm.nih.gov/Traces/sra/) under the accession numbers SRX495513 (Z5) and SRX495634 (M), while the unigene sequences data were also deposited in the Transcriptome Shotgun Assembly (TSA) project (PRJNA264349) at GenBank under the accession GBUF00000000. The version described in this paper is the first version, GBUF01000000. GO, eggNOG and KEGG classifications were further performed for the unigenes. In GO analysis, 6,699 unigenes of the total 19,186 (34.92%) were categorised into three major functional groups ‘biological process’, ‘cellular component’ and ‘molecular function’ and 4,718 unigenes (62.37%) were assigned to more than one GO terms. With the eggNOG analysis, a total of 14,218 unigenes of the total 19,186 (74.11%) were categorised into 25 groups (Fig. 1). Except genes with ‘Functional unknown (20.17%)’, 3,179 (15.50%) unigenes were assigned to ‘Signal transduction mechanisms’, 2433 (11.09%) unigenes to ‘General function prediction only’ and 1733 (7.90%) unigenes to ‘Transcription’ (Fig. 1). According to KEGG enrichment analysis, a total of 12,720 unigenes were consequently classified into specific pathways. Among them, 3,135 unigenes (24.64%) were included in the basic metabolism processes and the most important processes were relevant to ‘Carbohydrate metabolism (729 unigenes)’, ‘Amino Acid Metabolism (477 unigenes)’, ‘Lipid Metabolism (377 unigenes)’ and ‘Energy Metabolism (354 unigenes)’ (Fig. 2). These data greatly enrich transcriptomic resource of E. sinensis at different developmental stages and were expected to facilitate further molecular study on crabs.

3.3. Differentially expressed genes (DEGs) between the two developmental stages

3.2. Sequences assembly and annotation Totally, 68,574,292 and 63,433,542 raw reads from Z5 and M of E. sinensis were generated. After filtering, 62,654,562 and 58,339,580 clean reads remained with 6.11 Gb and 5.70 Gb of data for Z5 and M, respectively. All clean reads generated 276,571 contigs and assembly of the contigs resulted in 133,001 transcripts with an average length of 949 bp, which is higher than the transcript number in the previous study for the first zoeae (100, 252; Cui et al., 2013). Characteristics of the data sets were summarised in Table 2. After the BlastX searches against the NR database, 19,186 unigenes were obtained with annotations (E-value b 1E-05). The range of the unigene length was from Table 2 Summary information of transcriptomic sequences of the fifth zoeae (Z5) and megalopae (M) of E. sinensis. Parameters

Z5

M

Raw reads (n) Q20 percentage (%) GC percentage (%) Clean reads (n) Clean bases (Gb) Average length of clean reads (bp)

68,574,292 97.09 51.33 62,654,562 6.11 97.50

63,433,542 97.35 51.64 58,339,580 5.70 97.65

Assembly statistics Total contigs (n) Average contig length (max) (bp) N50 of contigs (bp) Total transcripts (n) Average transcript length (max) (bp) N50 of transcripts (bp) Annotated unigenes (n) Average unigene length (max) (bp) N50 of unigenes (bp)

276,571 339.82 (27,999) 545 133,001 949 (18,964) 2005 19,186 1,697 (18,964) 2571

The expression patterns of DEGs between these two stages were then identified. There were 652 unigenes from the total 19,186 unigenes (3.40% of the unigenes) with significant difference (p b 0.05) between Z5 and M stages, including 524 up-regulated genes at Z5 stage and 128 up-regulated genes at M stage (Supplementary file 2). The distribution of the significant changes was illustrated in a volcano plot (Fig. 3). To understand the functional distribution of these DEGs, all the DEGs were further analysed by GO and KEGG enrichment, respectively. Based on GO enrichment analysis, these 652 DEGs were also categorised into three major functional groups, ‘biological process’, ‘cellular component’ and ‘molecular function’. In the ‘biological process’, 44 subcategories were included (Fig. 4). Among these groups, ‘oxidoreductase activity (11%)’, ‘small molecule metabolic process (7.89%)’, ‘cellular nitrogen compound metabolic process (7.89%)’, ‘peptidase activity (6.81%)’, ‘biosynthetic process (6.81%)’, ‘catabolic process (5.73%)’, ‘transport (4.66%)’, ‘carbohydrate metabolic process’, ‘anatomical structure development (4.30%)’ and ‘signal transduction (4.23%)’ were the top ten processes with larger proportions of DEGs between the two stages (Fig. 4). KEGG enrichment analysis revealed that 305 DEGs were involved in 35 predicated metabolic pathways (Supplementary file 3). The most important metabolic pathways that might play roles in the larval development of E. sinensis from Z5 to M included ‘amino acid metabolism (20.27%)’, ‘carbohydrate metabolism (20.24%)’, ‘lipid metabolism (16.21%)’, ‘energy metabolism (8.78%)’, nucleotide metabolism (8.78%) and metabolism of cofactors and vitamins (8.78%) (Supplementary file 3). The predicted pathways, together with the GO analysis, reveal that the ‘signal transduction’ and ‘metabolisms’ are the principally changed processes in the development. These results are useful for the following discovery of genes related to the predominant change from Z5 to M.

4

Y. Li et al. / Comparative Biochemistry and Physiology, Part D 13 (2015) 1–9

Fig. 1. eggNOG functional distribution of all unigenes in the transcriptomes of Eriocheir sinensis. (For interpretation of the references to colour in this figure, the reader is referred to the web version of this article.)

3.4. DEGs involved in the morphological and sensorium development of the larvae During the metamorphosis from Z5 to M, the morphology of the larvae changes dramatically. Capabilities of swimming, crawling, sensory and predation are greatly improved with the development of statocyst, pereiopods, stalked eyes and claws, respectively. In the ‘signal transduction’ pathways, the principal ones were the ‘MAPK signaling pathway’, ‘Wnt signaling pathway’ and ‘Hedgehog signaling pathway’(Supplementary file 3), which were revealed to be related to the periodic morphological changes of the prawn Macrobrachium

rosenbergii (Ventura et al., 2013). Except genes involved in the pathways, many other DEGs were also identified to possibly contribute to the morphological change. More than 30 cuticle protein genes were found expressed differently between the two stages. It demonstrates that during the metamorphosis, various copies of cuticle protein were involved in the metamorphosis as revealed in the crab Portunus pelagicus by moult cycle specific differential gene expression profiling (Kuballa et al., 2011). Other DEGs related to the formation or degradation of cytoskeleton and mussel were also identified. The seven up-regulated genes in Z5 were mucin-2 (MUC2), zinc metalloproteinase (MP), endocuticle structural glycoprotein SgAbd-2 (CUD2), lysin motif-containing protein (Lysin), actin,

Fig. 2. Functional distribution of all unigenes in the transcriptomes of Eriocheir sinensis based on KEGG analysis. (For interpretation of the references to colour in this figure, the reader is referred to the web version of this article.)

Y. Li et al. / Comparative Biochemistry and Physiology, Part D 13 (2015) 1–9

5

Fig. 3. Volcano plot of differentially expressed genes (DEGs) from the transcriptomes of the fifth zoeae (Z5) and megalopae (MA) in Eriocheir sinensis. For each unigene, the ratio of expression levels (Z5 V.S. M) was plotted against the −log error rate. The horizontal line indicates the significance threshold (p b 0.05), and the vertical lines indicate the two fold change threshold. Non-differentially expressed genes are shown with orange dots, and DEGs are shown with blue dots. (For interpretation of the references to colour in this figure, the reader is referred to the web version of this article.)

helix–loop–helix protein Delilah (HLH) and myosin IIIA (Myosin), while the two highly expressed genes in M were articulin p60 (articulin p60) and envelope protein (Envelope) (Table 3). In the development process, two

hormone-related genes, ecdysteroid-regulated 16 kDa protein (ESR16) and hormone receptor 4-like isoform 1 (HR4-1), were up-regulated in Z5, and one hormone transport gene was detected to be up-regulated

Fig. 4. GO distributions of the differentially expressed genes involved in biological processes from the transcriptomes of the fifth zoeae and megalopae in Eriocheir sinensis. (For interpretation of the references to colour in this figure, the reader is referred to the web version of this article.)

6

Y. Li et al. / Comparative Biochemistry and Physiology, Part D 13 (2015) 1–9

Table 3 Differentially expressed genes related to the morphological and sensory perception changes from the fifth zoeae (Z5) to megalopae (M) of Eriocheir sinensis. Gene name

Protein name

Matched organism

Fold change (M:Z5)⁎

Morphological change Mucin comp44879_c0_seq1 MP comp35516_c0_seq1 CUD2 comp34350_c0_seq1 ESR16 comp42970_c0_seq1 Lysin comp45903_c1_seq1 Actin comp24445_c0_seq2 HR4-1 comp48176_c0_seq2 HLH comp47099_c0_seq4 Myosin comp49034_c0_seq1 AUX comp49129_c1_seq1 Articulin comp38708_c0_seq10 Envelope comp53892_c0_seq1

Unigene ID

Mucin-2 Zinc metalloproteinase Endocuticle structural glycoprotein SgAbd-2 Ecdysteroid-regulated 16 kDa protein precursor Lysin motif-containing protein Actin Hormone receptor 4-like isoform 1 Helix–loop–helix protein delilah Myosin IIIA Auxin transport protein Articulin p60 Envelope protein

Ceratitis capitata Saccoglossus kowalevskii Ceratitis capitata Bombyx mori Procambarus clarkii Mayetiola destructor Acyrthosiphon pisum Tribolium castaneum Danaus plexippus Tribolium castaneum Rimicaris exoculata Cotesia congregata

0 0.008 0.082 0.092 0.103 0.123 0.145 0.157 0.181 10.424 19.944 Inf

Sensory perception RLBP1 BCMO1 POE SDR16C5 CHRNA10 ACE-1 UGT8

Retinaldehyde-binding protein 1 beta-Carotene 15,15'-monooxygenase 1 Protein purity of essence Epidermal retinal dehydrogenase 2 Nicotinic acetylcholine receptor subunit alpha 10 Acetylcholinesterase 1 2-Hydroxyacylsphingosine 1-beta-galactosyltransferase

Camponotus floridanus Salmo salar Acromyrmex echinatior Crassostrea gigas Pandalopsis japonica Liposcelis entomophila Pteropus alecto

0.037 0.109 5.780 6.917 8.297 8.517 11.831

comp39285_c0_seq1 comp49781_c0_seq1 comp49129_c0_seq1 comp46769_c0_seq5 comp47167_c7_seq1 comp50739_c0_seq1 comp47096_c0_seq4

⁎ Numbers in bold means genes expression is up-regulated in M and the rest up-regulated in Z5; ‘Inf’ means genes only expressed in M. The expression of genes in red color is validated by qRT-PCR.

in M, which was auxin transport protein (AUX) (Table 3). Among them, ESR16 and HR4-1 can coordinates growth and maturation during larval development and play central roles in the genetic cascades triggered by the steroid hormone ecdysone at the onset of metamorphosis. It is therefore evident that the molting and metamorphosis in the mitten crab is genetically controlled through endocrine systems that mediate gene expression as reviewed (McWilliam and Phillips, 2007). As the development, the expression of genes related to sensory capability and nervous system displayed significant differences between Z5 and M. In Z5, two up-regulated genes involved in visual development were identified, which were retinaldehyde-binding protein 1 (RLBP1) and beta-carotene 15,15'-monooxygenase 1 (BCMO1) (Table 3). It is noted that RLBP1 was a functional component of the visual cycle, and knockout of retinaldehyde binding protein in mice result in delayed dark adaptation (Saari et al., 2001). In M, five up-regulated genes were identified to be related to the sensory perception, including protein purity of essence (POE), epidermal retinal dehydrogenase 2 (SDR16C5), nicotinic acetylcholine receptor subunit alpha 10 (CHRNA10), acetylcholinesterase 1 (ACE1) and 2-hydroxyacylsphingosine 1-beta-galactosyltransferase (UGT8) (Table 3). Among them, POE and SDR16C5 are related to visual sense; CHRNA10 is relevant to auditory stimulus, while ACE and UGT8 are involved in the development of nervous system. The expression patterns of the DEGs demonstrate that the sensory and nervous systems become more mature with the development, such as the formation of large stalked eyes in M, which might enhance the activity capability of the M.

3.5. DEGs related to food digestion and nutrition metabolism Although previous research has revealed many genes involved in the food digestion and nutrition metabolism in crabs (Jiang et al., 2009; Wang et al., 2014), very limited molecular information is available for understanding these processes in the larvae. It is known that the larvae at M stage were more euryphagous than Z5. M can not only filter feeding plankton but also prey on the large zooplankon, such as cladocerans and copepods. Moreover, larvae at M stage are more ferocious and agile, which can easily capture Z5. Therefore, the genes related to food digestion and nutrition metabolism are expected to be expressed in different patterns in the larvae at the two stages. There were 62 DEGs related to digestive absorption and substance metabolism being identified from the 652 DEGs, with 35 genes up-regulated in Z5 and 27 upregulated in M. They were divided into the following three categories (Tables 4, 5, 6). 3.5.1. Carbohydrate metabolism Although some researches show that growth was not affected by the carbohydrate content (Bages and Sloane, 1981), it is believed that carbohydrate plays an important role in launching the utilisation of protein and lipid for energy production (Rutledge, 1999; Holme et al., 2009), while carbohydrase is the key for the carbohydrate digestion from diet (Perera et al., 2008; Jiang et al., 2009). It was proven that the expression of carbohydrase in P. trituberculatus was different at various

Table 4 Differentially expressed genes associated with carbohydrate metabolism between the fifth zoeae (Z5) and the megalopae (M) of Eriocheir sinensis. Gene name

Unigene ID

Protein name

Matched organism

Fold change (M:Z5)⁎

GDH GDH SnEG54 CHST12 TPI GXYLT1 CHST11 B3GNT1 GLYS C2orf18

comp41799_c1_seq2 comp42263_c1_seq2 comp12693_c0_seq1 comp39922_c0_seq1 comp37480_c0_seq1 comp48491_c0_seq1 comp37429_c1_seq1 comp47217_c1_seq1 comp44526_c0_seq1 comp48222_c0_seq2

Glucose dehydrogenase Glucose dehydrogenase Cellulase Chondroitin 4-sulfotransferase Triosephosphate isomerase Glucoside xylosyltransferase 1 Carbohydrate sulfotransferase 11 N-acetyllactosaminide beta-1,3-N-acetylglucosaminyltransferase Glycogen synthase Transmembrane protein C2orf18

Megachile rotundata Aedes aegypti Mesocentrotus nudus Aedes aegypti Penaeus monodon Strongylocentrotus purpuratus Nasonia vitripennis Megachile rotundata Danaus plexippus Ixodes scapularis

0.004 0.026 0.065 0.079 0.081 0.087 0.129 0.130 5.964 7.002

⁎ Numbers in bold means genes expression is up-regulated in M and the rest up-regulated in Z5; ‘Inf’ means genes only expressed in M.

Y. Li et al. / Comparative Biochemistry and Physiology, Part D 13 (2015) 1–9

7

Table 5 Differentially expressed genes associated with lipid metabolism between the fifth zoeae (Z5) and megalopae (M) of Eriocheir sinensis. Gene name

Unigene ID

Protein name

Matched organism

Fold change (M:Z5)⁎

LIP3 PNLIP PLA2 PLA2 LIP1 APOD ELOVL FAS ACOT2 MID1ip SCD9 ACADS ACOX3 AGMO AGPAT

comp41807_c0_seq1 comp48856_c0_seq10 comp47494_c4_seq7 comp43539_c0_seq2 comp40112_c0_seq1 comp49330_c0_seq2 comp45991_c1_seq1 comp50119_c0_seq1 comp31966_c0_seq1 comp46211_c0_seq1 comp29189_c0_seq1 comp42005_c0_seq1 comp47909_c0_seq1 comp31174_c0_seq1 comp46750_c1_seq1

Lipase 3 Triacylglycerol lipase, pancreatic Phospholipase A2 Phospholipase A2, major isoenzyme Triacylglycerol lipase Apolipoprotein D Elongation of very long chain fatty acids protein Fatty acid synthase Acyl-coenzyme A thioesterase 2, mitochondrial Mid1-interacting protein Acyl-CoA delta-9 desaturase Short-chain specific acyl-CoA dehydrogenase, mitochondrial Peroxisomal acyl-coenzyme A oxidase 3 Alkylglycerol monooxygenase 1-Acyl-sn-glycerol-3-phosphate acyltransferase alpha

Nasonia vitripennis Danaus plexippus Harpegnathos saltator Trichechus manatus latirostris Litopenaeus vannamei Harpegnathos saltator Marsupenaeus japonicus Litopenaeus vannamei Cricetulus griseus Ixodes scapularis Eriocheir sinensis Apis florea Takifugu rubripes Gallus gallus Acromyrmex echinatior

0.013 0.014 0.020 0.020 0.021 0.046 0.104 5.967 6.764 8.482 9.626 10.194 11.238 46.939 51.731

⁎ Numbers in bold means genes expression is up-regulated in M and the rest up-regulated in Z5; ‘Inf’ means genes only expressed in M. The expression of genes in red color is validated by qRT-PCR.

developmental stages (Wang et al., 2014). In our study, 10 DEGs were identified to be related to carbohydrate digestion and metabolism. Eight genes were up-regulated in Z5, including two carbohydrate sulfotransferases (CHSTs), two glucose dehydrogenases (GDHs) one triosephosphate isomerase (TPI) and three other carbohydrate digestion-related genes (Table 4), while only two genes related to carbohydrate metabolism were up-regulated in M, which were glycogen synthase (GLYS) and transmembrane protein C2orf18 (C2orf18) (Table 4). It indicates that the larvae at Z5 stage might require more carbohydrate than M, which might be used to compensate for the energy consumed in metamorphosis from Z5 to M. As reported (Rutledge, 1999; Johnston, 2003), carbohydrates are more important for early development stage of crustaceans than late development stage. It is therefore obvious that the ability to utilise carbohydrate in the diet varies not only among species, but also at different larval developmental stages in the crabs. 3.5.2. Lipid metabolism In the lipid metabolism, 15 unigenes were detected to be differentially expressed between Z5 and M, with seven DEGs up-regulated in Z5 and eight up-regulated in M (Table 5). Among them, more lipid degradation-related genes were up-regulated in Z5, e.g., lipase 3

(LIP3), triacylglycerol lipase (PNLIP) and phospholipase A2 (PLA2), while more fatty acid synthases were found highly expressed in M, such as fatty acid synthase (FAS), Mid1-interacting protein (MID1ip), acyl-CoA delta-9 desaturase (SCD9) and peroxisomal acyl-coenzyme A oxidase 3 (ACOX3). Decomposing of the fatty acid is the primary source yielding energy, and therefore our results suggest that more energy is possibly required in Z5, while in M, more energy is reserved in the form of fatty acids. It might be explained that during the metamorphosis, more energy is essential to be released for Z5 to complete the transition into M. 3.5.3. Protein metabolism Protein is the key group of essential nutrients required by all animals for growth. The knowledge of the protein metabolism is therefore important for studying the development of a species. In this study, it was shown that in the 20 DEGs related to the protein digestion and metabolism, 11 genes were up-regulated in Z5 (e.g., trypsin, PRSS; aminopeptidase N, ANPEP; chymotrypsin-like elastase family member 2A, CELA2A) and nine genes were up-regulated in M (e.g., allantoicase 1, ALLB1; glycine dehydrogenase, GLDC; peroxisomal sarcosine oxidase, PIPOX) (Table 6). Among these, both trypsin and chymotrypsin belong to serine proteases and it was noted that one other serine protease 27 (SP) gene

Table 6 Differentially expressed genes associated with protein metabolism between the fifth zoeae (Z5) and megalopae (M) of Eriocheir sinensis. Gene name

Unigene ID

Protein name

Matched organism

Fold change (M:Z5)⁎

PRSS8 PRSS PHGDH GCP SP ANPEP DPF-6 CELA2A PAM WAP2 ALR ALLB1 UROC1 HGDO TRY GLDC PIPOX HPD BHMT SP22D

comp41136_c0_seq2 comp44886_c0_seq1 comp47822_c0_seq4 comp42772_c0_seq1 comp42784_c0_seq1 comp46290_c0_seq1 comp29390_c0_seq1 comp46794_c1_seq1 comp39954_c1_seq2 comp35782_c0_seq1 comp48055_c0_seq1 comp46108_c0_seq3 comp46860_c1_seq1 comp43696_c0_seq1 comp45079_c3_seq3 comp42180_c0_seq1 comp49434_c0_seq1 comp44606_c0_seq3 comp42600_c0_seq3 comp35785_c1_seq1

Prostasin Trypsin d-3-Phosphoglycerate dehydrogenase Plasma glutamate carboxypeptidase Serine protease 27 Minopeptidase N Dipeptidyl peptidase family member 6 Chymotrypsin-like elastase family member 2A Peptidylglycine alpha-hydroxylating monooxygenase WAP2 Alanine racemase Allantoinase 1 Urocanate hydratase 1 Homogentisate 1,2-dioxygenase Tyrosine/tryptophan monooxygenase Glycine dehydrogenase (decarboxylating) Peroxisomal sarcosine oxidase 4-Hydroxyphenylpyruvate dioxygenase Betaine-homocysteine methyltransferase Serine protease 22D

Takifugu rubripes Pediculus humanus corporis Strongylocentrotus purpuratus Apis florea Oreochromis niloticus Harpegnathos saltator Crassostrea gigas Apis florea Procambarus clarkii Eriocheir sinensis Marsupenaeus japonicus Crassostrea gigas Xenopus (Silurana) tropicalis Pediculus humanus corporis Pediculus humanus corporis Xenopus laevis Maylandia zebra Aedes aegypti Saccoglossus kowalevskii Anopheles gambiae

0.003 0.009 0.014 0.020 0.032 0.037 0.045 0.075 0.137 0.145 0.163 5.749 6.039 6.383 7.726 8.588 9.109 15.281 17.219 18.139

⁎ Numbers in bold means genes expression is up-regulated in M and the rest up-regulated in Z5; ‘Inf’ means genes only expressed in M. The expression of genes in red color is validated by qRT-PCR.

8

Y. Li et al. / Comparative Biochemistry and Physiology, Part D 13 (2015) 1–9

was also found expressed highly in Z5 (Table 6), suggesting their important roles in the protein digestion, which were also documented in the study for P. trituberculatus (Wang et al., 2014). Accordingly, our results indicate that the protein metabolism process is active in Z5 and more protein might be digested at the stage. In the previous study, most diets formulated for crustacean larvae contain 30–50% crude protein (McConaugha, 1985; Le Moullac and Van Wormhoudt, 1994). According to our results, less protein is required in M than Z5, and it was similar as previous researches. In those studies, it was reported that the appropriate protein requirements at Z5 and M were 45% and 40%, respectively, and the survival and metamorphic rate of E. sinensis larvae changed significantly when they were fed different protein content (Pan et al., 2005; Sui et al., 2007). Therefore, our study also provides useful molecular information for the aquaculture breeding.

3.6. qRT-PCR verification The DEGs were selected to verify the results of the RNA-Seq analysis by qRT-PCR, using different RNA from those used for RNA-Seq. All the tested genes showed significant differentially expression between Z5 and M (Fig. 5). In qRT-PCR, Mucin, RLBP1, PLA2 and PRSS showed low abundance in M, while AGMO and SP22D showed up-regulation in M. Even though most qRT-PCR results (Fig. 5) indicated smaller differences (except AGMO) compared with the RNA-Seq analysis (in red color; Tables 3, 5 and 6), there was a consistent expression tendency between the two results. This verifies that the accuracy of the RNA-Seq. In conclusion, this study is the first investigation into the transcriptomic change between the fifth zoeae (Z5) and the megalopae (M) of E. sinensis, which are the key larvae developmental stages in the

Fig. 5. qRT-PCR verification of differentially expressed unigenes. Z5 means the fifth zoeal sample and F means the megalopal sample. *Statistically significant difference is detected (p b 0.05). (For interpretation of the references to colour in this figure, the reader is referred to the web version of this article.)

Y. Li et al. / Comparative Biochemistry and Physiology, Part D 13 (2015) 1–9

crab. The differentially expressed genes (DEGs) between the two stages were identified, especially genes relevant to morphological change, sensory ability, neuroendocrine system and food digestion as well as nutrition metabolism. However, many transcripts identified as differentially regulated are functional unknown, which require further study. Nevertheless, our study provides useful information for the future molecular genetic study on the issues. Supplementary data to this article can be found online at http://dx. doi.org/10.1016/j.cbd.2014.10.002. Acknowledgements This work was supported by grants from the Chinese National ‘863’ Project (no. 2012AA10A409) and the National Natural Science Foundation of China (31101924 and 41276165) to Dr. Zhaoxia Cui. References Anders, S., Huber, W., 2010. Differential expression analysis for sequence count data. Genome Biol. 11, R106. Bages, M., Sloane, L., 1981. Effects of dietary protein and starch levels on growth and survival of Penaeus monodon (Fabricius) postlarvae. Aquaculture 25, 117–128. Benjamini, Y., Hochberg, Y., 1995. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Stat. Soc. B 57, 289–300. Biesiot, P.M., Capuzzo, J.M., 1990. Changes in digestive enzyme activities during early development of the American lobster Homarus americanus Milne Edwards. J. Exp. Mar. Biol. Ecol. 136, 107–122. Cheng, Y., Wu, X., Yang, X., Hines, A.H., 2008. Current trends in hatchery techniques and stock enhancement for Chinese mitten crab, Eriocheir japonica sinensis. Rev. Fish. Sci. 16, 377–384. Conesa, A., Gotz, S., Garcia-Gomez, J.M., Terol, J., Talon, M., Robles, M., 2005. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics 21, 3674–3676. Cox, M.P., Peterson, D.A., Biggs, P.J., 2010. SolexaQA: at-a-glance quality assessment of Illumina second-generation sequencing data. BMC Bioinforma. 11, 485. Cui, Z.X., Li, X.H., Liu, Y., Song, C.W., Hui, M., Shi, G.H., Luo, D.L., Li, Y.D., 2013. Transcriptome profiling analysis on whole bodies of microbial challenged Eriocheir sinensis larvae for immune gene identification and SNP development. PLoS One 8, e82156. Gonçalves, A.T., Farlora, R., Gallardo-Escárate, C., 2014. Transcriptome survey of the lipid metabolic pathways involved in energy production and ecdysteroid synthesis in the salmon louse Caligus rogercresseyi (Crustacea: Copepoda). Comp. Biochem. Physiol. B 176, 9–17. Gosselin, L.A., Qian, P.Y., 1997. Juvenile mortality in benthic marine invertebrates. Oceanogr. Lit. Rev. 44, 730–731. Grabherr, M.G., Haas, B.J., Yassour, M., Levin, J.Z., Thompson, D.A., Amit, I., Adiconis, X., Fan, L., Raychowdhury, R., Zeng, Q.D., Chen, Z.H., Mauceli, E., Hacohen, N., Gnirke, A., Rhind, N., di Palma, F., Birren, B.W., Nusbaum, C., Lindblad-Toh, K., Friedman, N., Regev, A., 2011. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat. Biotechnol. 29, 644–652. He, L., Wang, Q., Jin, X.K., Wang, Y., Chen, L.L., Liu, L.H., Wang, Y., 2012. Transcriptome profiling of testis during sexual maturation stages in Eriocheir sinensis using Illumina sequencing. PLoS One 7, e33735. He, L., Jiang, H., Cao, D.D., Liu, L.H., Hu, S.N., Wang, Q., 2013. Comparative transcriptome analysis of the accessory sex gland and testis from the Chinese mitten crab (Eriocheir sinensis). PLoS One 8. Herborg, L.M., Rushton, S., Clare, A., Bentley, M., 2003. Spread of the Chinese mitten crab (Eriocheir sinensis H. Milne Edwards) in Continental Europe: analysis of a historical data set. Hydrobiologia 503, 21–28. Holme, M.-H., Zeng, C., Southgate, P.C., 2009. A review of recent progress toward development of a formulated microbound diet for mud crab, Scylla serrata, larvae and their nutritional requirements. Aquaculture 286, 164–175. Hunt, H.L., Scheibling, R.E., 1997. Role of early post-settlement mortality in recruitment of benthic marine invertebrates. Mar. Ecol. Prog. Ser. 155, 269–301. Jiang, H., Yin, Y.X., Zhang, X.W., Hu, S.N., Wang, Q., 2009. Chasing relationships between nutrition and reproduction: A comparative transcriptome analysis of hepatopancreas and testis from Eriocheir sinensis. Comp. Biochem. Physiol. D 4, 227–234. Johnston, D., 2003. Ontogenetic changes in digestive enzyme activity of the spiny lobster, Jasus edwardsii (Decapoda; Palinuridae). Mar. Biol. 143, 1071–1082. Kanehisa, M., Goto, S., Kawashima, S., Okuno, Y., Hattori, M., 2004. The KEGG resource for deciphering the genome. Nucleic Acids Res. 32, D277–D280.

9

Kuballa, A.V., Holton, T.A., Paterson, B., Elizur, A., 2011. Moult cycle specific differential gene expression profiling of the crab Portunus pelagicus. BMC Genomics 12, 147. Laufer, H., Biggers, W.J., 2001. Unifying concepts learned from methyl farnesoate for invertebrate reproduction and post-embryonic development. Am. Zool. 41, 442–457. Laufer, H., Borst, D., Baker, F., Reuter, C., Tsai, L., Schooley, D., Carrasco, C., Sinkus, M., 1987. Identification of a juvenile hormone-like compound in a crustacean. Science 235, 202–205. Le Moullac, G., Van Wormhoudt, A., 1994. Adaptation of digestive enzymes to dietary protein, carbohydrate and fibre levels and influence of protein and carbohydrate quality in Penaeus uannamei larvae (Crustacea, Decapoda). Aquat. Living Resour. 7, 203–210. Lemos, D., Hernández-Cortés, M., Navarrete, A., Garcia-Carreno, F., Phan, V., 1999. Ontogenetic variation in digestive proteinase activity of larvae and postlarvae of the pink shrimp Farfantepenaeus paulensis (Crustacea: Decapoda: Penaeidae). Mar. Biol. 135, 653–662. Li, X.H., Cui, Z.X., Liu, Y., Song, C.W., Shi, G.H., 2013. Transcriptome analysis and discovery of genes involved in immune pathways from hepatopancreas of microbial challenged mitten crab Eriocheir sinensis. PLoS One 8, e68233. Lovett, D.L., Felder, D.L., 1990. Ontogenetic change in digestive enzyme activity of larval and postlarval white shrimp Penaeus setiferus (Crustacea, Decapoda, Penaeidae). Biol. Bull. 178, 144–159. McConaugha, J.R., 1985. Nutrition and larval growth. In: Wenner, A.M. (Ed.), CrustaceanLarval Growth 2. A.A. Balkema, Rotterdam, pp. 127–154. McWilliam, P.S., Phillips, B.F., 2007. Spiny lobster development: mechanisms inducing metamorphosis to the puerulus: a review. Rev. Fish Biol. Fish. 17, 615–632. Montú, M., Anger, K., De Bakker, C., 1996. Larval development of the Chinese mitten crab Eriocheir sinensis H. Milne-Edwards (Decapoda: Grapsidae) reared in the laboratory. Helgoländer Meeresun. 50, 223–252. Mortazavi, A., Williams, B.A., Mccue, K., Schaeffer, L., Wold, B., 2008. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat. Methods 5, 621–628. Osman, R.W., Whitlatch, R.B., 2004. The control of the development of a marine benthic community by predation on recruits. J. Exp. Mar. Biol. Ecol. 311, 117–145. Palmer, M.A., Allan, J.D., Butman, C.A., 1996. Dispersal as a regional process affecting the local dynamics of marine and stream benthic invertebrates. Trends Ecol. Evol. 11, 322–326. Pan, L.Q., Xiao, G.Q., Zhang, H.X., Luan, Z.H., 2005. Effects of different dietary protein content on growth and protease activity of Eriocheir sinensis larvae. Aquaculture 246, 313–319. Perera, E., Moyano, F., Díaz, M., Perdomo-Morales, R., Montero-Alejo, V., Alonso, E., Carrillo, O., Galich, G., 2008. Polymorphism and partial characterization of digestive enzymes in the spiny lobster Panulirus argus. Comp. Biochem. Physiol. B 150, 247–254. Powell, S., Szklarczyk, D., Trachana, K., Roth, A., Kuhn, M., Muller, J., Arnold, R., Rattei, T., Letunic, I., Doerks, T., Jensen, L.J., von Mering, C., Bork, P., 2012. eggNOG v3.0: orthologous groups covering 1133 organisms at 41 different taxonomic ranges. Nucleic Acids Res. 40, D284–D289. Roughgarden, J., Gaines, S., Possingham, H., 1988. Recruitment dynamics in complex life cycles. Science 241, 1460–1466. Rutledge, S., 1999. Characterisation of cellulase activity and implications for diet formulation in the mud crab, Scylla serrata. Honours Thesis, School of Natural Resource Sciences. Queensland University of Technology, Brisbane, Qld, Australia (65 pp.). Saari, J.C., Nawrot, M., Kennedy, B.N., Garwin, G.G., Hurley, J.B., Huang, J., Possin, D.E., Crabb, J.W., 2001. Visual cycle impairment in cellular retinaldehyde binding protein (CRALBP) knockout mice results in delayed dark adaptation. Neuron 29, 739–748. Shen, H.S., Hu, Y.C., Zhou, X., 2014. Sex-lethal gene of the Chinese mitten crab Eriocheir sinensis: cDNA cloning, induction by eyestalk ablation, and expression of two splice variants in males and females. Dev. Genes Evol. 224, 97–105. Sui, L., Wille, M., Cheng, Y., Sorgeloos, P., 2007. The effect of dietary n-3 HUFA levels and DHA/EPA ratios on growth, survival and osmotic stress tolerance of Chinese mitten crab Eriocheir sinensis larvae. Aquaculture 273, 139–150. Sui, L., Wille, M., Cheng, Y., Wu, X., Sorgeloos, P., 2011. Larviculture techniques of Chinese mitten crab Eriocheir sinensis. Aquaculture 315, 16–19. Ventura, T., Manor, R., Aflalo, E.D., Chalifa-Caspi, V., Weil, S., Sharabi, O., Sagi, A., 2013. Post-Embryonic transcriptomes of the prawn Macrobrachium rosenbergii: Multigenic succession through metamorphosis. PLoS One 8, e55322. Wahle, R.A., 2003. Revealing stock-recruitment relationships in lobsters and crabs: is experimental ecology the key? Fish. Res. 65, 3–32. Wang, W., Wu, X., Liu, Z., Zheng, H., Cheng, Y., 2014. Insights into hepatopancreatic functions for nutrition metabolism and ovarian development in the crab Portunus trituberculatus: gene discovery in the comparative transcriptome of different hepatopancreas stages. PLoS One 9, e84921. Zhang, W., Wan, H.L., Jiang, H., Zhao, Y.L., Zhang, X.W., Hu, S.N., Wang, Q., 2011. A transcriptome analysis of mitten crab testes (Eriocheir sinensis). Genet. Mol. Biol. 34, 136–141.