- Email: [email protected]
Comparative analysis of transcriptomes from different coloration of Chinese mitten crab Eriocheir sinensis
Fish and Shellfish Immunology 98 (2020) 515–521
Contents lists available at ScienceDirect
Fish and Shellfish Immunology journal homepage: www.elsevier.com/locate/fsi
Full length article
Comparative analysis of transcriptomes from different coloration of Chinese mitten crab Eriocheir sinensis
T
Yingying Zhaoa, Xiaochen Zhua, Zhibin Hana, Yazhao Zhanga, Tengfei Donga, Yingdong Lia, Jing Donga, Hua Weia, Xiaodong Lia,b,∗ a b
Key Laboratory of Zoonosis of Liaoning Province, College of Animal Science and Veterinary Medicine, Shenyang Agricultural University, Shenyang, 110866, China Panjin Guanghe Crab Industry Co.Ltd., Panjin, 124000, China
A R T I C LE I N FO
A B S T R A C T
Keywords: Eriocheir sinensis Comparative transcriptome Different expressed genes Immune-related gene Color-related gene
Chinese mitten crab Eriocheir sinensis is probably the most important freshwater cultured crab in China. A tiny minority of brownish-orange individuals have been discovered in the long period of artificial breeding history of E. sinensiss. Those mutants are usually accompanied with slow growth rate, low molting frequency and poor survival rate, which may be the results of growth defects and immunodeficiency. To better understand the relationship between body color determination and the immune system as well as the related genes expression in E. sinensiss, we performed the whole-body transcriptome analysis in different color of first stage zoea (ZI) larvae using next-generation sequencing (NGS) technology. We randomly assembled 175.40 and 177.52million clean reads from the wild and mutant ZIs, respectively. Finally, we identified 7153 differentially expressed genes (DEGs) (p < 0.05), with 5194 up-regulated and 1959 down-regulated. A total of 13 KEGG pathways related to immune system were detected among 248 pathways. Except the first whole-body RNA sequencing of colorspecific transcriptomes for E. sinensis, this study will offer a better understanding of the underlying molecular mechanisms of interaction between color determination and the immune system.
1. Introduction Shell color, which can be involved in communication, reproduction, thermoregulation, and self-protection from predators, plays an important role in the life history of crustaceans [1]. Except species-specific color difference, intraspecific color deviation, which varies with body size, sex [2] or locality [3], also has been found in some crustaceans. Being different from momentary body color shift as a result of dispersion or aggregation of melanophore pigment [4,5], phenotypic color differences may directly affect its market price and ornamental value. For example, darker-colored black tiger prawns (Penaeus monodon) often has higher commodity price in Australian market [6]. Therefore, body color is often used as a phenotypic marker for breeding of some economic crustaceans, such as Geothelphusa dehaani [7], Penaeus monodon [6] and Fenneropenaeus merguiensis [8], and have attracted many scientists to explore the underlying mechanisms of coloration and color pattern formation [9]. Hitherto, color-related gene expression has been well studied in some aquatic animals, such as red tilapia [10], sea cucumbers [11,12], and pacific oyster [13], etc. As for crustacean, in addition to the environmental background and nutrients that influence
their shell color [14–16], Apo-crustacyanin has been considered as an essential color-related gene in prawn and lobsters genetically [17–19]. Chinese mitten crab, Eriocheir sinensis, is probably the most important freshwater cultured crab in China [20]. Although adult E. sinensiss are normally greenish-brown in color, a tiny minority of brownish-orange individuals have been discovered in the long period of artificial breeding history (Fig. 1). Those individuals exhibit apparent red from heart beating stage of zoea to juvenile, and gradually turn into brownish-orange after sexual maturity. It is predictable that the brownish-orange individuals have higher market value than those with normal color due to their bright color. Unfortunately, significantly lower specific growth rate and molting frequency as well as poor survival rate have been observed in brownishorange individuals in the comparison with normal colored counterparts [21], which may be the results of growth defects and immunodeficiency. However, those mutants provide an excellent experimental material to find out the formation mechanism of E. sinensis body color, and to explore the relationship between crustacean body color and growth, molting, immune system and so on. Next-generation sequencing (NGS) technologies have been widely
∗ Corresponding author. Key Laboratory of Zoonosis of Liaoning Province, College of Animal Science and Veterinary Medicine, Shenyang Agricultural University, Shenyang, 110866, China. E-mail addresses: [email protected] (Y. Zhao), [email protected] (X. Li).
https://doi.org/10.1016/j.fsi.2020.01.051 Received 8 December 2019; Received in revised form 21 January 2020; Accepted 24 January 2020 Available online 27 January 2020 1050-4648/ © 2020 Elsevier Ltd. All rights reserved.
Fish and Shellfish Immunology 98 (2020) 515–521
Y. Zhao, et al.
Fig. 1. Two body color types of Eriocheir sinensis. (A, mutant type; B, wild type). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)
used to generate a mass of gene expression information and transcript sequences for non-model organisms without genome data [22,23]. So far, many color-related genes have been reported in fish and mollusks through comparative transcriptomes [10–13]. However, comparative transcriptome for different colored crustacean remain unexplored. This study was designed to spot new candidate genes involved in color variation, and to elaborate the different expression patterns between different color types of E. sinensis. Besides, immunity related DEGs were determined and analyzed. Our outcomes may be helpful for the enrichment of genetic information and further study on underlying molecular mechanisms of color variation in E. sinensis.
org/) for possible protein domain. In addition, all assembled transcripts were searched against the NCBI non-redundant Protein database (NR), String Protein database, Swiss-Prot protein database, Cluster of Orthologous Groups database (COG) as well as Kyoto Encyclopedia of Genes and Genomes pathway database (KEGG), KEGG Orthologue database (KO), and Gene Ontology (GO) using BlastX with an Evalue = 1E-05. The best hit species were identified based on NR annotation, and GO annotations were determined using Blast2GO to obtain the functional classification of the unigenes [25].
2. Materials and methods
The unigene expression was measured in accordance with the fragments per kilobase of exon region per million mapped reads (FPKM) method using RSEM [26]. Differentially expressed genes (DEGs) were identified by the edgeR program [27], and the selection criteria was false discovery rate (FDR) < 0.05 and |log2(fold change)|≥ 1. Eventually, GO and KEGG analyses of DEGs were performed by Goatools (https://github.com/tanghaibao/GOatools) [28], and KOBAS (http://kobas.cbi.pku.edu.cn/home.do) [29], respectively.
2.4. Differential expression analysis
2.1. Sample preparation Healthy wild (greenish-brown) and mutant (brownish-orange) berried females of E. sinensis were collected from Panjin Guanghe Crab Industry Co. Ltd. in Liaoning Province, China. They were reared under the same conditions before sampling. Three lumps of samples (approximate 2 mL in a lump), which were in first stage zoea (ZI), were collected from one greenish-brown and one brownish-orange female, respectively. Each lump of ZI larvae was separately stored in liquid nitrogen for subsequent RNA extraction.
2.5. Experimental verification by quantitative real-time PCR To validate the reliability of RNA-seq results, new wild and mutant E. sinensis ZI larvae were collected from other berried females. Ten putative genes (Astacin, Lipocalin, CaM 1, CTSL, HSP70, LSH-R, MBP, MEP1A, CHH1, PKA C) were chosen for verification of gene expression between different body color types of E. sinensis by quantitative realtime PCR (qRT-PCR) analysis. Primer Premier 3.0 software (http:// bioinfo.ut.ee/primer3-0.4.0/) was used for designing PCR primers (Table S1). Reverse-transcriptions of cDNA templates were carried out using the HiScript Q RT reagent Kit. The SYBR Green RT-PCR assay was performed by the ABI PRISM 7500 Sequence Detection System (Applied Biosystems). The relative expression levels were normalized by reference gene (β-actin) through 2−ΔΔCt method. At last, the data were analyzed by ABI 7500 SDS software (Applied Biosystems). For every gene, qRT-PCR reactions for each sample were repeated three times. SPSS 22.0 (IBM Corp. Released 2013. IBM SPSS Statistics for Windows, Version 22.0. Armonk, NY: IBM Corp.) was used for one-way analysis of variance and the results were shown as mean ± S.D.
2.2. RNA extraction and sequencing For sequencing, total RNA extraction of six different samples was performed separately according to the protocol associated of TRIzol® Reagent Kit (Invitrogen, USA). The RNA integrality and quality were examined by agarose gel electrophoresis and Agilent 2100 Bioanalyzer (Agilent, USA). The mRNA was purified using Dynabeads (Invitrogen, USA), and the final concentration of mRNA was assessed using NanoDrop2000 (Thermo Scientific, USA). Six libraries (3 for greenish-brown and 3 for brownish-orange) of E. sinensis were obtained. The cDNA libraries were constructed by using the Illumina Truseq™ RNA sample prep Kit (Illumina, USA) following the manufacturer's instruction. The libraries were sequenced using Illumina HiSeq4000 instrument by Shanghai Majorbio Bio-pharm Technology Corporation, China. All clean Illumina sequencing data have been uploaded to SRA database.
3. Results and discussion 2.3. De novo assembly and gene annotation 3.1. Transcriptome sequencing and assembly After data filtering, poor quality reads such as adaptors, low quality sequences (< Q20), and the reads less than 20 bp in length were removed. The remaining reads were de novo assembled by software Trinity (http://trinityrnaseq.sourceforge.net/) [24]. All assembled sequences were searched against the Pfam database (http://pfam.sanger.ac.uk/) using HMMER3 program (http://hmmer.
Six cDNA libraries were generated, and a total of 175.40 and 177.52million clean reads for the wild and mutant ZI larvae were obtained, respectively (Table 1, SRA: PRJNA595668). The sequencing quality was high, with Q20 ratio > 96% and Q30 ratio > 90% of each sample. And the average GC content was 51.96% (Table 1). 516
Fish and Shellfish Immunology 98 (2020) 515–521
Y. Zhao, et al.
Table 1 Summary of transcriptome sequencing and assembly of E. sinensis. Parameters
Reads Total raw reads (Mb) Total Clean reads (Mb) Total Clean bases (Gb) Clean reads Q20 (%) Clean reads Q30 (%) GC percentage (%)
Assembly statistics Total Number Average length (bp) GC percentage (%) N50 of transcripts (bp)
Wild groups
Mutational groups
ES1
ES2
ES3
ESR1
ESR2
ESR3
64.39 61.63 8.80 96.67 90.50 52.82
55.60 53.47 7.67 96.70 90.54 53.12
62.42 60.26 8.67 96.85 90.89 52.44
59.06 56.51 8.08 96.71 90.61 51.10
59.87 57.59 8.24 96.85 90.98 51.04
65.55 63.42 9.11 97.00 91.30 51.21
transcripts 82360 1307 49.45 2249
unigenes 63387 1065 48.85 1758
After de novo assembly using Trinity, 82360 transcripts and 63387 unigenes were assembled with an average length of 1307bp and 1065bp, respectively. The GC content of these assembled transcripts was 49.45%, and N50 length was 2249 bp (Table 1). These values were higher than previous study of E. sinensis transcriptome-analysis results [30]. The most of the transcripts were less than 800 bp (51.18%), while 19.13% transcripts were over 2000 bp. These outcomes indicated that the sequences data of present study were high-quality and reliable due to its high quality, and more information on relevant genes were expected to be collected.
Fig. 2. Venn diagram of unigene numbers annotated in five databases.
single-organism process (3750, 40.6%) and regulation of biological process (1549, 16.8%) subcategories in biological process category. Besides, 66 unigenes were enriched in immune system process (0.7%) (Fig. S1). 6538 unigenes were categorized into 25 categories in KOG database (Fig. S2, Table S2). Among them, the largest number of unigenes were assigned into the categories of translation, ribosomal structure and biogenesis (915, 14.0%), followed by signal transduction mechanisms (774, 11.8%) and posttranslational modification, protein turnover and chaperones (611, 9.3%), successively. On the contrary, the smallest numbers of unigenes were assigned into extracellular structures (20, 0.3%), nuclear structure (10, 0.2%) and cell motility (5, 0.1%). In addition, 24 unigenes (0.4%) were assigned to defense mechanisms. After KEGG pathway analysis, 17132 unigenes were classified into 372 known pathways and the category with the most unigenes is Ribosome (1263, 7.3%), followed by Endocytosis (526, 3.1%), Huntington's disease (495, 2.9%), Protein processing in endoplasmic reticulum (461, 2.7%), and RNA transport (450, 2.6%). The 372 pathways were grouped into five biochemical pathways, including metabolism (A), genetic information processing (B), environmental information processing (C), cellular processes (D) and organismal systems (E) (Fig. 4). “Global and overview maps” (2671, 15.6%), was the most abundant subcategories within the category of metabolism. As for genetic information processing, environmental information processing and cellular processes, the most common subcategory was “Translation” (2382, 13.9%), “Signal transduction” (1942, 11.3%) and “Cell growth and death” (1026, 6.0%), respectively. “Endocrine system” (1221, 7.1%) and “immune system” (983, 5.7%) were the two most abundant subcategories in organismal systems (Table S3).
3.2. Functional annotation and classification Total of 63387 unigenes were annotated with mean length of 1065 bp, and N50 length of 1758 bp were further assembled. Similar with the transcripts, most unigenes were less than 800 bp (58.48%), and many transcripts were more than 2000 bp in length (13.66%). All of the unigenes were annotated in at least one database, and 3373 genes (5.32%) of them were annotated in all five databases. The annotation success rate of unigenes was 22015 in NR (34.73%), 17132 in KEEG (27.03%), 20919 in SwissProt (33.00%), 17316 in Pfam (27.32%), and 6538 in String (10.31%), respectively (Table 2, Fig. 2). The species distribution revealed that E. sinensis has the highest number of hits to the Zootermopsis nevadensis (1086, 8.2%), followed by Daphnia pulex (1290, 5.8%), Tribolium castaneum (745, 3.4%), Pediculus humanus (500, 2.3%), Nasonia vitripennis (371, 1.7%) and other species (Fig. 3). 3.3. Functional classification, enrichment, and pathway analysis 9243 of the 12236 predicted transcripts were successfully assigned to GO annotation within the three primary GO categories (molecular functional, cellular component, and biological process) and 59 sub-categories based on the NR annotation. Most unigenes were assigned to metabolic process (5273, 57.0%), cellular process (5133, 55.5%),
3.4. Identification and analysis of DEGs
Table 2 Statistics of unigene functional annotation. Annotation in Database
Unigene No.
Percentage (%)
NR KEGG Swissprot Pfam String All annotated transcripts
22015 17132 20919 17316 6538 3373
34.73 27.03 33.00 27.32 10.31 5.32
After the filtration by the fold changes, 7153 DEGs were identified between the wild and mutant ZI larvae, with 5194 of them were upregulated and 1959 were down-regulated. Statistical significance as well as magnitude of change was visualized by volcano plot (Fig. S3). The degree in change of up-regulated DEGs was apparently higher than that of down-regulated DEGs. GO function and pathway enrichment analysis divided all DEGs into three categories, biological process, cellular component and molecular 517
Fish and Shellfish Immunology 98 (2020) 515–521
Y. Zhao, et al.
Fig. 3. Comparison of E. sinensis transcriptomic sequences with the known sequences in database.
membrane-bounded organelle”, and “non-membrane-bounded organelle” (p < 0.001). KEGG pathway analysis assigned a total of 3539 DEGs into 248 known pathways, which belonged to five core categories, cellular process, metabolism, environmental information processing, genetic information processing and organismal systems. “Transport and catabolism” (223, 45.9%) in the category of cellular process, “Carbohydrate metabolism” (240, 22.4%) in the category of metabolism, “Signal transduction” (394, 87.6%) in the category of environmental
function. The highest percentage of DEGs involved in GO terms was “metabolic process” of the biological processes (613, 8.6%), “cellular component” (607, 8.5%), in the category of cellular component, and “molecular function” (930, 13.0%) in the category of molecular functions. Finally, the top ten significantly changed categories were as followed, “ribonucleoprotein complex biogenesis”, “structural constituent of ribosome”, “cellular macromolecule biosynthetic process”, “ribosome biogenesis”, “structural molecule activity”, “ribosome”, “translation”, “macromolecule biosynthetic process”, “intracellular non-
Fig. 4. KEGG functional classification of assembled unigenes. 518
Fish and Shellfish Immunology 98 (2020) 515–521
Y. Zhao, et al.
Table 3 Immune-related pathways in KEGG and Go Terms. Database KEGG ko04612 ko04640 ko04670 ko04622 ko04611 ko04621 ko04610 ko04666 ko04660 ko04062 ko04620 ko04662 ko04650 GO GO:0002697 GO:0002698 GO:0002683
Pathway
Up No.
Down No.
Total No.
p value
Antigen processing and presentation Hematopoietic cell lineage Leukocyte transendothelial migration RIG-I-like receptor signaling pathway Platelet activation NOD-like receptor signaling pathway Complement and coagulation cascades Fc gamma R-mediated phagocytosis T cell receptor signaling pathway Chemokine signaling pathway Toll-like receptor signaling pathway B cell receptor signaling pathway Natural killer cell mediated cytotoxicity
33 2 10 1 13 6 1 14 3 11 0 2 1
6 4 8 3 8 2 1 14 2 11 3 0 1
39 6 18 4 21 8 2 28 5 22 2 2 2
0.00501 0.28483 0.46669 0.51429 0.52441 0.63878 0.70776 0.92131 0.96407 0.99027 0.99087 0.99449 0.99821
regulation of immune effector process negative regulation of immune effector process negative regulation of immune system process
2 1 1
1 1 1
3 2 2
0.00526 0.0127 0.0353
information processing, the “Translation” (386, 66.3%) subcategory of “genetic information processing” category, and “Cardiac muscle contraction” (53) of organismal systems category was the most common subcategories, respectively. Moreover, the most abundant pathways included “Ribosome” (293, 23.2%), “Oxidative phosphorylation” (89, 22.6%), “Lysosome” (76, 21.3%) and “Carbon metabolism” (67, 15.0%) and “Protein processing in endoplasmic reticulum” (66, 14.3%). The top five significantly distinct categories were “Ribosome”, “Oxidative phosphorylation”, “Lysosome”, “Cardiac muscle contraction” and “Proteasome” (p < 0.001).
ATPase subunit alpha (ATP), glutamate receptor 1 (AMPAR), glutamate receptor (NMDAR), cAMP-dependent protein kinase catalytic subunit (PKA), Ras-related protein Rap-1 (Rap1), Ras-like GTP-binding protein RHO (Rho), serine/threonine-protein phosphatase PP1 catalytic subunit (PP1), Aconitate hydratase (ACO), p21-activated kinase 1 (PAK1), phosphatidylinositol phospholipase C epsilon (PLCε) and Ras homolog gene family member A (RhoA) in the cAMP signaling pathway, as well as low density lipoprotein receptor-related protein 5/6 (LRP5/6), casein kinase II subunit alpha (CK2), glycogen synthase kinase 3 beta (GSK3β), PKA, RhoA, serine/threonine-protein phosphatase 2B catalytic subunit (CaN) and C-terminal binding protein (CtBP) in the Wnt signaling pathway. Most of those DEGs exhibited significantly up-regulated expression. It is well known that melanin has various functions in different organisms, such as protecting cell from UV light damage for human, improving avian feathers’ resistance to abrasion, etc. [42]. Beyond that, for crustacean, melanin plays a significant role in immune response by involving in wound melanization, encapsulation of foreign particles, clot formation, and so on [43,44]. Although the underlying mechanisms are still unclear, those DEGs may make a significant contribution to the shell color variation of E. sinensis and the immune deficiencies of brownish-orange E. sinensis. Moreover, GO analysis identified 7 DEGs enriched in the “Immune system process” including “regulation of immune effector process”, “negative regulation of immune effector process” and “negative regulation of immune system process” (p < 0.05) (Table 3). Specifically, one of lipocalins, Maltose-binding periplasmic protein (MBP), was enriched in all the three processes. As transport proteins, lipocalin family takes part in many biological functions, including retinol transport, pheromone transport, prostaglandin synthesis, invertebrate camouflage, olfaction, immune response, etc. [45]. In the present study, 34 DEGs were related to the lipocalins, and 13 of them were significantly different, including the Apolipoprotein D-like, crustacyanin-like lipocalin, putative cellular retinoic acid/retinol binding protein (CRABP-II), MBP, Fatty acid-binding protein and Apolipoprotein D (ApoD). Among them, the Apolipoprotein D-like, crustacyaninlike lipocalin, and MBP were all down-regulated whereas CRABP-II, Fatty acid-binding protein and ApoD were up-regulated in mutant E. sinesis. This may illustrate that up/down-regulation of lipocalin family proteins were related to body color of E. sinesis.
3.5. Analysis of DEGs of immune-related genes Finally, 13 KEGG pathways related to immune system were detected among 248 pathways mapped to the KEGG database. The numbers of DEGs in the 13 pathways were listed in Table 3. Our KEGG analysis showed that the up/down-regulated DEGs were present in 13 “immune system” pathways, and enrichment in most pathways was not significant except pathway “Antigen processing and presentation” (p < 0.05) (Table 3). 33 up-regulated and 6 down-regulated DEGs were enriched on ten unigenes, including HSP70 (heat shock protein 70), HSP90 (heat shock protein 90), CTSB/L/S (cathepsin B/L/S) and so on (Fig. S4). HSP can protect crustacean from various internal or external stresses, and it is evidently involved in innate immune response [31–34]. CTSB was reported to participate in the regulation of developmental processes in some arthropods [35]. In addition, as an important participant in the MHC class I pathway, most of the proteasome subunits are significantly up-regulated in mutant type (Fig. S5). The proper expression of proteasome plays an important role in cellular homeostasis, and proteasome dysregulation have been associated with many human diseases, such as cancer and Parkinson's disease [36–39]. Similarly, abnormal proteasome expression was also found in mud crab and silkworm in their response to viral or fungal infection [40,41]. In the present study, combined with the overexpression of heat shock protein and proteasome, it may infer that the mutant individuals are experiencing some kind of stress (probably caused by genetic defect). In this study, with the exception of pathway “Antigen processing and presentation”, some significantly differential DEGs have been located in melanogenesis related pathways, such as cyclic adenosine monophosphate/protein kinase A (cAMP/PKA) signaling pathway, Wnt signaling pathway, mitogen-activated protein kinases (MAPK) signaling pathway, and so on. For example, we found DEGs including adrenergic receptor beta-1(GPCR), adenylate cyclase 1 (AC), cyclic nucleotide gated channel alpha 1 (CNGC), Calmodulin (CaM), Plasma membrane calcium-transporting ATPase (PMCA), sodium/potassium-transporting
3.6. Validation of ten DEGs by qRT-PCR The expression of ten chosen DEGs were verified on three individual biological replicates through qRT-PCR. Among ten chosen genes, four of them were down-regulated in mutant type (Astacin, Lipocalin, MBP 519
Fish and Shellfish Immunology 98 (2020) 515–521
Y. Zhao, et al.
Appendix A. Supplementary data Supplementary data to this article can be found online at https:// doi.org/10.1016/j.fsi.2020.01.051. References [1] N.G. Ertl, A. Elizur, P. Brooks, A.V. Kuballa, T.A. Anderson, W.R. Knibb, Molecular characterisation of colour formation in the prawn Fenneropenaeus merguiensis, PloS One 8 (2013) e56920, , https://doi.org/10.1371/journal.pone.0056920. [2] Y. Murakami, K. Wada, Inter-populational variations in body color related to growth stage and sex in Gaetice depressus (De Haan, 1835) (Decapoda, Brachyura, Varunidae), Crustaceana 88 (2015) 113–126, https://doi.org/10.1163/1568540300003391. [3] A.S. Freire, M.A.A. Pinheiro, H. Karam-Silva, M.M. Teschima, Biology of Grapsus grapsus (Linnaeus, 1758) (Brachyura, Grapsidae) in the Saint Peter and Saint Paul Archipelago, Equatorial Atlantic ocean, Helgol. Mar. Res. 65 (2011) 263–273, https://doi.org/10.1007/s10152-010-0220-5. [4] M.Z. Darnell, Ecological physiology of the circadian pigmentation rhythm in the fiddler crab Uca panacea, J. Exp. Mar. Biol. Ecol. 426–427 (2012) 39–47, https:// doi.org/10.1016/j.jembe.2012.05.014. [5] S.M. Kronstadt, M.Z. Darnell, P. Munguia, Background and temperature effects on Uca panacea color change, Mar. Biol. 160 (2013) 1373–1381, https://doi.org/10. 1007/s00227-013-2189-5. [6] R.K. Tume, A.L. Sikes, S. Tabrett, D.M. Smith, Effect of background colour on the distribution of astaxanthin in black tiger prawn (Penaeus monodon): effective method for improvement of cooked colour, Aquaculture 296 (2009) 129–135, https://doi.org/10.1016/j.aquaculture.2009.08.006. [7] M. Ikeda, T. Suzuki, Y. Fujio, Genetic differentiation among populations of Japanese freshwater crab, Geotheiphusa dehaani (White), with reference to the body color variation, Benthos Res. 53 (1998) 47–52, https://doi.org/10.5179/benthos1996. 53.1_47. [8] N.H. Nguyen, J. Quinn, D. Powell, A. Elizur, N.P. Thoa, J. Nocillado, R. Lamont, C. Remilton, W. Knibb, Heritability for body colour and its genetic association with morphometric traits in Banana shrimp (Fenneropenaeus merguiensis), BMC Genet. 15 (2014) 1–12, https://doi.org/10.1186/s12863-014-0132-5. [9] J.J. Nordlund, R.E. Boissy, V.J. Hearing, R.A. King, W.S. Oetting, J.P. Ortonne, The Pigmentary System: Physiology and Pathophysiology, second ed., Wiley–Blackwell, 2007 1–1229. [10] W. Zhu, L. Wang, Z. Dong, X. Chen, F. Song, N. Liu, H. Yang, J. Fu, Comparative transcriptome analysis identifies candidate genes related to skin color differentiation in red tilapia, Sci. Rep. 6 (2016) 1–12, https://doi.org/10.1038/srep31347. [11] D. Ma, H. Yang, L. Sun, D. Xu, Comparative analysis of transcriptomes from albino and control sea cucumbers, Apostichopus japonicus, Acta Oceanol. Sin. 33 (2014) 55–61, https://doi.org/10.1007/s13131-014-0464-z. [12] L. Xing, L. Sun, S. Liu, X. Li, L. Zhang, H. Yang, Transcriptome analysis provides insights into the mechanism of albinism during different pigmentation stages of the albino sea cucumber Apostichopus japonicus, Aquaculture 486 (2018) 148–160, https://doi.org/10.1016/j.aquaculture.2017.12.016. [13] D. Feng, Q. Li, H. Yu, X. Zhao, L. Kong, Comparative transcriptome analysis of the pacific oyster Crassostrea gigas characterized by shell colors: identification of genetic bases potentially involved in pigmentation, PloS One 10 (2015) 1–17, https://doi. org/10.1371/journal.pone.0145257. [14] A.T. Palma, R.S. Steneck, Does variable coloration in juvenile marine crabs reduce risk of visual predation ? Ecology 82 (2001) 2961–2967, https://doi.org/10.2307/ 2679974. [15] M.F. Tlusty, A. Metzler, S. Huckabone, S. Suanda, S. Guerrier, Morphological colour change in the American lobster (Homarus americanus) in response to background colour and UV light, N. Z. J. Mar. Freshw. Res. 43 (2009) 247–255, https://doi.org/ 10.1080/00288330909509998. [16] N.M. Wade, S. Cheers, N. Bourne, S. Irvin, D. Blyth, B.D. Glencross, Dietary astaxanthin levels affect colour, growth, carotenoid digestibility and the accumulation of specific carotenoid esters in the Giant Tiger Shrimp, Penaeus monodon, Aquacult. Res. 48 (2017) 395–406, https://doi.org/10.1111/are.12888. [17] M. Wang, X. Zhu, J. Yang, Z. Dai, K. Mahmood, F. Yang, W. Yang, Prawn lipocalin: characteristics and expressional pattern in subepidermal adipose tissue during reproductive molting cycle, Comp. Biochem. Physiol. B Biochem. Mol. Biol. 147 (2007) 222–229, https://doi.org/10.1016/j.cbpb.2007.01.003. [18] N.M. Wade, A. Tollenaere, M.R. Hall, B.M. Degnan, Evolution of a novel carotenoidbinding protein responsible for crustacean shell color, Mol. Biol. Evol. 26 (2009) 1851–1864, https://doi.org/10.1093/molbev/msp092. [19] F. Yang, M. Wang, Y. Ma, W. Ma, W. Yang, Prawn lipocalin: characterization of a color shift induced by gene knockdown and ligand binding assay, J. Exp. Zool. Part A Ecol. Genet. Physiol. 315 A (2011) 562–571, https://doi.org/10.1002/jez.706. [20] D. Chen, M. Zhang, S. Shrestha, Compositional characteristics and nutritional quality of Chinese mitten crab (Eriocheir sinensis), Food Chem. 103 (2007) 1343–1349, https://doi.org/10.1016/j.foodchem.2006.10.047. [21] Y. Si, X. Li, Y. Jiang, Y. Wang, X. Liu, Y. Zheng, J. Suo, N. Sun, Culture pattern of juvenile Chinese mitten handed crab Eriocheir sinensis with red color shell, Journal of Dalian Ocean University 32 (2017) 139–144, https://doi.org/10.16535/j.cnki. dlhyxb.2017.02.003. [22] M.L. Metzker, Sequencing technologies the next generation, Nat. Rev. Genet. 11 (2010) 31–46, https://doi.org/10.1038/nrg2626. [23] M. Garber, M.G. Grabherr, M. Guttman, C. Trapnell, Computational methods for
Fig. 5. qRT-PCR verification of differentially expressed genes between the wild and mutant E. sinensis. *Significant level (p < 0.05). ** (p < 0.01),*** (p < 0.001).
and CHH1), whereas other six were up-regulated (CaM 1, CTSL, HSP70, LSH-R, MEP1A, PKA C) (Fig. 5). The results of qRT-PCR were consistent with expression pattern of RNA-Seq analyses, which support the accuracy and reliability of DEGs given by RNA-Seq. 4. Conclusion In conclusion, we constructed and sequenced the color-specific transcriptomes ZI larvae of E. sinensis. Illumina paired-end RNA sequencing technology was used to sequence. In total, More than 350 million clean reads were obtained, which then assembled into 82360 transcripts and 63387 unigenes. The unigenes were annotated by six databases, NR, Pfam, KOG, Swiss-prot, KEGG, and GO, and we obtained comprehensive gene function information. We identified 7153 significantly DEGs including 5194 significantly up-regulated and 1959 down-regulated DEGs. KEGG and GO analysis indicated those DEGs in pathway “Antigen processing and presentation”, melanogenesis related pathway, as well as those of lipocalin family protein may associated with shell color alteration and immunodeficiency of mutant E. sinensis. Overall, this study provides valuable information on the molecular mechanisms underlying interactions between color determination and the immune system in E. sinensis, but also facilitates further investigations of functional genomics for this species and other closely related species. Author contributions Y Zhao and X Li conceived and designed the experiments. Z Han, Y Zhang, T Dong, H Wei, J Dong and Y Li performed the experiments. Y Zhao and X Zhu analyzed the data and wrote the paper. X Li contributed the reagents and materials. Acknowledgments This work was supported by grants from the Cultivation Plan for Youth Agricultural Science and Technology Innovative Talents of Liaoning Province (No. 2015044), Talent Introduction Program of Shenyang Agricultural University (No. 880416005), Guide Project for Key Research and Development Program of Liaoning province (No. 201802120), and Science and Technology Plan Project of Liaoning Province (No. 2019JH2/10200006). 520
Fish and Shellfish Immunology 98 (2020) 515–521
Y. Zhao, et al.
[24]
[25]
[26]
[27]
[28]
[29]
[30]
[31]
[32]
[33]
237–246, https://doi.org/10.1016/j.fsi.2016.11.049. [34] X. Chu, T. Yang, Y. Liu, L. Hong, T. Jiao, X. Meng, D. Zhang, J. Wang, B. Tang, C. Zhou, Q. Liu, W. Zhang, W. He, Transcriptome analysis of differential expressed genes in hepatopancreas of Procambarus clarkii challenged with peptidoglycan, Fish Shellfish Immunol. 86 (2019) 311–318, https://doi.org/10.1016/j.fsi.2018.11.048. [35] H. Wang, C. Fu, Q. Zeng, F. Li, H. Wang, J. Sun, Comparative transcriptome analysis reveals related regulatory mechanisms of androgenic gland in Eriocheir sinensis, BioMed Res. Int. 1 (2017) 1–12, https://doi.org/10.1155/2017/4956216. [36] A. Kumatori, K. Tanaka, N. Inamura, S. Sone, T. Ogura, T. Matsumoto, T. Tachikawa, S. Shin, A. Ichihara, Abnormally high expression of proteasomes in human leukemic cells, Proc. Natl. Acad. Sci. U. S. A 87 (1990) 7071–7075, https:// doi.org/10.1073/pnas.87.18.7071. [37] L. Chen, K. Madura, Increased proteasome activity, ubiquitin-conjugating enzymes, and eEF1A translation factor detected in breast cancer tissue, Canc. Res. 65 (2005) 5599–5606, https://doi.org/10.1158/0008-5472.CAN-05-0201. [38] Q. Zheng, T. Huang, L. Zhang, Y. Zhou, H. Luo, H. Xu, X. Wang, Dysregulation of ubiquitin-proteasome system in neurodegenerative diseases, Front. Aging Neurosci. 8 (2016) 1–10, https://doi.org/10.3389/fnagi.2016.00303. [39] R. Motosugi, S. Murata, Dynamic regulation of proteasome expression, Front. Mol. Biosci. 6 (2019) 4–11, https://doi.org/10.3389/fmolb.2019.00030. [40] C. Hou, G. Qin, T. Liu, T. Geng, K. Gao, Z. Pan, H. Qian, X. Guo, Transcriptome analysis of silkworm, Bombyx mori, during early response to Beauveria bassiana challenges, PloS One 9 (2014) 1–13, https://doi.org/10.1371/journal.pone. 0091189. [41] S. Liu, G. Chen, H. Xu, W. Zou, W. Yan, Q. Wang, H. Deng, H. Zhang, G. Yu, J. He, S. Weng, Transcriptome analysis of mud crab (Scylla paramamosain) gills in response to Mud crab reovirus (MCRV), Fish Shellfish Immunol. 60 (2017) 545–553, https://doi.org/10.1016/j.fsi.2016.07.033. [42] F. Solano, Melanins: skin pigments and much more—types, structural models, biological functions, and formation routes, New J. Sci. 2014 (2014) 1–28, https:// doi.org/10.1155/2014/498276. [43] H. Bilandžija, M. Laslo, M.L. Porter, D.W. Fong, Melanization in response to wounding is ancestral in arthropods and conserved in albino cave species, Sci. Rep. 7 (2017) 1–11, https://doi.org/10.1038/s41598-017-17471-2. [44] L. Vazquez, J. Alpuche, G. Maldonado, C. Agundis, A. Pereyra-Morales, E. Zenteno, Immunity mechanisms in crustaceans, Innate Immun. 15 (2009) 179–188, https:// doi.org/10.1177/1753425909102876. [45] D.R. Flower, The lipocalin protein family: structure and function, Biochem. J. 318 (1996) 1–14, https://doi.org/10.1042/bj3180001.
transcriptome annotation and quantification using RNA-seq, Nat. Methods 8 (2011) 469–477, https://doi.org/10.1038/nmeth.1613. M.G. Grabherr, B.J. Haas, M. Yassour, J.Z. Levin, D.A. Thompson, I. Amit, X. Adiconis, L. Fan, R. Raychowdhury, Q. Zeng, Z. Chen, E. Mauceli, N. Hacohen, A. Gnirke, N. Rhind, F. Di Palma, B.W. Birren, C. Nusbaum, K. Lindblad-Toh, N. Friedman, A. Regev, Full-length transcriptome assembly from RNA-Seq data without a reference genome, Nat. Biotechnol. 29 (2011) 644–652, https://doi.org/ 10.1038/nbt.1883. A. Conesa, S. Götz, J.M. García-Gómez, J. Terol, M. Talón, M. Robles, Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research, Bioinformatics 21 (2005) 3674–3676, https://doi.org/10.1093/ bioinformatics/bti610. B. Li, C.N. Dewey, RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome, BMC Bioinf. 12 (2011) 323–338, https://doi. org/10.1093/10.1186/1471-2105-12-323. M.D. Robinson, D.J. McCarthy, G.K. Smyth, edgeR: a Bioconductor package for differential expression analysis of digital gene expression data, Bioinformatics 26 (2009) 139–140, https://doi.org/10.1093/bioinformatics/btp616. D.V. Klopfenstein, L. Zhang, B.S. Pedersen, F. Ramírez, A.W. Vesztrocy, A. Naldi, C.J. Mungall, J.M. Yunes, O. Botvinnik, M. Weigel, W. Dampier, C. Dessimoz, P. Flick, H. Tang, GOATOOLS: a Python library for Gene Ontology analyses, Sci. Rep. 8 (2018) 1–17, https://doi.org/10.1038/s41598-018-28948-z. C. Xie, X. Mao, J. Huang, Y. Ding, J. Wu, S. Dong, L. Kong, G. Gao, C. Li, L. Wei, KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases, Nucleic Acids Res. 39 (2011) 316–322, https://doi.org/10.1093/nar/ gkr483. Y. Liu, M. Hui, Z. Cui, D. Luo, C. Song, Y. Li, L. Liu, Comparative transcriptome analysis reveals sex-biased gene expression in juvenile Chinese mitten crab Eriocheir Sinensis, PloS One 10 (2015) 1–21, https://doi.org/10.1371/journal.pone.0133068. L. Cerenius, P. Jiravanichpaisal, H. Liu, I. Söderhäll, Crustacean Immunity, Invertebrate Immunity, Springer, Boston, MA, 2010, pp. 239–259, https://doi.org/ 10.1007/978-1-4419-8059-5_13. M.K. Chaurasia, F. Nizam, G. Ravichandran, M.V. Arasu, N.A. Al-Dhabi, A. Arshad, P. Elumalai, J. Arockiaraj, Molecular importance of prawn large heat shock proteins 60, 70 and 90, Fish Shellfish Immunol. 48 (2016) 228–238, https://doi.org/10. 1016/j.fsi.2015.11.034. W. Junprung, P. Supungul, A. Tassanakajon, HSP70 and HSP90 are involved in shrimp Penaeus vannamei tolerance to AHPND-causing strain of Vibrio parahaemolyticus after non-lethal heat shock, Fish Shellfish Immunol. 60 (2017)
521
Contents lists available at ScienceDirect
Fish and Shellfish Immunology journal homepage: www.elsevier.com/locate/fsi
Full length article
Comparative analysis of transcriptomes from different coloration of Chinese mitten crab Eriocheir sinensis
T
Yingying Zhaoa, Xiaochen Zhua, Zhibin Hana, Yazhao Zhanga, Tengfei Donga, Yingdong Lia, Jing Donga, Hua Weia, Xiaodong Lia,b,∗ a b
Key Laboratory of Zoonosis of Liaoning Province, College of Animal Science and Veterinary Medicine, Shenyang Agricultural University, Shenyang, 110866, China Panjin Guanghe Crab Industry Co.Ltd., Panjin, 124000, China
A R T I C LE I N FO
A B S T R A C T
Keywords: Eriocheir sinensis Comparative transcriptome Different expressed genes Immune-related gene Color-related gene
Chinese mitten crab Eriocheir sinensis is probably the most important freshwater cultured crab in China. A tiny minority of brownish-orange individuals have been discovered in the long period of artificial breeding history of E. sinensiss. Those mutants are usually accompanied with slow growth rate, low molting frequency and poor survival rate, which may be the results of growth defects and immunodeficiency. To better understand the relationship between body color determination and the immune system as well as the related genes expression in E. sinensiss, we performed the whole-body transcriptome analysis in different color of first stage zoea (ZI) larvae using next-generation sequencing (NGS) technology. We randomly assembled 175.40 and 177.52million clean reads from the wild and mutant ZIs, respectively. Finally, we identified 7153 differentially expressed genes (DEGs) (p < 0.05), with 5194 up-regulated and 1959 down-regulated. A total of 13 KEGG pathways related to immune system were detected among 248 pathways. Except the first whole-body RNA sequencing of colorspecific transcriptomes for E. sinensis, this study will offer a better understanding of the underlying molecular mechanisms of interaction between color determination and the immune system.
1. Introduction Shell color, which can be involved in communication, reproduction, thermoregulation, and self-protection from predators, plays an important role in the life history of crustaceans [1]. Except species-specific color difference, intraspecific color deviation, which varies with body size, sex [2] or locality [3], also has been found in some crustaceans. Being different from momentary body color shift as a result of dispersion or aggregation of melanophore pigment [4,5], phenotypic color differences may directly affect its market price and ornamental value. For example, darker-colored black tiger prawns (Penaeus monodon) often has higher commodity price in Australian market [6]. Therefore, body color is often used as a phenotypic marker for breeding of some economic crustaceans, such as Geothelphusa dehaani [7], Penaeus monodon [6] and Fenneropenaeus merguiensis [8], and have attracted many scientists to explore the underlying mechanisms of coloration and color pattern formation [9]. Hitherto, color-related gene expression has been well studied in some aquatic animals, such as red tilapia [10], sea cucumbers [11,12], and pacific oyster [13], etc. As for crustacean, in addition to the environmental background and nutrients that influence
their shell color [14–16], Apo-crustacyanin has been considered as an essential color-related gene in prawn and lobsters genetically [17–19]. Chinese mitten crab, Eriocheir sinensis, is probably the most important freshwater cultured crab in China [20]. Although adult E. sinensiss are normally greenish-brown in color, a tiny minority of brownish-orange individuals have been discovered in the long period of artificial breeding history (Fig. 1). Those individuals exhibit apparent red from heart beating stage of zoea to juvenile, and gradually turn into brownish-orange after sexual maturity. It is predictable that the brownish-orange individuals have higher market value than those with normal color due to their bright color. Unfortunately, significantly lower specific growth rate and molting frequency as well as poor survival rate have been observed in brownishorange individuals in the comparison with normal colored counterparts [21], which may be the results of growth defects and immunodeficiency. However, those mutants provide an excellent experimental material to find out the formation mechanism of E. sinensis body color, and to explore the relationship between crustacean body color and growth, molting, immune system and so on. Next-generation sequencing (NGS) technologies have been widely
∗ Corresponding author. Key Laboratory of Zoonosis of Liaoning Province, College of Animal Science and Veterinary Medicine, Shenyang Agricultural University, Shenyang, 110866, China. E-mail addresses: [email protected] (Y. Zhao), [email protected] (X. Li).
https://doi.org/10.1016/j.fsi.2020.01.051 Received 8 December 2019; Received in revised form 21 January 2020; Accepted 24 January 2020 Available online 27 January 2020 1050-4648/ © 2020 Elsevier Ltd. All rights reserved.
Fish and Shellfish Immunology 98 (2020) 515–521
Y. Zhao, et al.
Fig. 1. Two body color types of Eriocheir sinensis. (A, mutant type; B, wild type). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)
used to generate a mass of gene expression information and transcript sequences for non-model organisms without genome data [22,23]. So far, many color-related genes have been reported in fish and mollusks through comparative transcriptomes [10–13]. However, comparative transcriptome for different colored crustacean remain unexplored. This study was designed to spot new candidate genes involved in color variation, and to elaborate the different expression patterns between different color types of E. sinensis. Besides, immunity related DEGs were determined and analyzed. Our outcomes may be helpful for the enrichment of genetic information and further study on underlying molecular mechanisms of color variation in E. sinensis.
org/) for possible protein domain. In addition, all assembled transcripts were searched against the NCBI non-redundant Protein database (NR), String Protein database, Swiss-Prot protein database, Cluster of Orthologous Groups database (COG) as well as Kyoto Encyclopedia of Genes and Genomes pathway database (KEGG), KEGG Orthologue database (KO), and Gene Ontology (GO) using BlastX with an Evalue = 1E-05. The best hit species were identified based on NR annotation, and GO annotations were determined using Blast2GO to obtain the functional classification of the unigenes [25].
2. Materials and methods
The unigene expression was measured in accordance with the fragments per kilobase of exon region per million mapped reads (FPKM) method using RSEM [26]. Differentially expressed genes (DEGs) were identified by the edgeR program [27], and the selection criteria was false discovery rate (FDR) < 0.05 and |log2(fold change)|≥ 1. Eventually, GO and KEGG analyses of DEGs were performed by Goatools (https://github.com/tanghaibao/GOatools) [28], and KOBAS (http://kobas.cbi.pku.edu.cn/home.do) [29], respectively.
2.4. Differential expression analysis
2.1. Sample preparation Healthy wild (greenish-brown) and mutant (brownish-orange) berried females of E. sinensis were collected from Panjin Guanghe Crab Industry Co. Ltd. in Liaoning Province, China. They were reared under the same conditions before sampling. Three lumps of samples (approximate 2 mL in a lump), which were in first stage zoea (ZI), were collected from one greenish-brown and one brownish-orange female, respectively. Each lump of ZI larvae was separately stored in liquid nitrogen for subsequent RNA extraction.
2.5. Experimental verification by quantitative real-time PCR To validate the reliability of RNA-seq results, new wild and mutant E. sinensis ZI larvae were collected from other berried females. Ten putative genes (Astacin, Lipocalin, CaM 1, CTSL, HSP70, LSH-R, MBP, MEP1A, CHH1, PKA C) were chosen for verification of gene expression between different body color types of E. sinensis by quantitative realtime PCR (qRT-PCR) analysis. Primer Premier 3.0 software (http:// bioinfo.ut.ee/primer3-0.4.0/) was used for designing PCR primers (Table S1). Reverse-transcriptions of cDNA templates were carried out using the HiScript Q RT reagent Kit. The SYBR Green RT-PCR assay was performed by the ABI PRISM 7500 Sequence Detection System (Applied Biosystems). The relative expression levels were normalized by reference gene (β-actin) through 2−ΔΔCt method. At last, the data were analyzed by ABI 7500 SDS software (Applied Biosystems). For every gene, qRT-PCR reactions for each sample were repeated three times. SPSS 22.0 (IBM Corp. Released 2013. IBM SPSS Statistics for Windows, Version 22.0. Armonk, NY: IBM Corp.) was used for one-way analysis of variance and the results were shown as mean ± S.D.
2.2. RNA extraction and sequencing For sequencing, total RNA extraction of six different samples was performed separately according to the protocol associated of TRIzol® Reagent Kit (Invitrogen, USA). The RNA integrality and quality were examined by agarose gel electrophoresis and Agilent 2100 Bioanalyzer (Agilent, USA). The mRNA was purified using Dynabeads (Invitrogen, USA), and the final concentration of mRNA was assessed using NanoDrop2000 (Thermo Scientific, USA). Six libraries (3 for greenish-brown and 3 for brownish-orange) of E. sinensis were obtained. The cDNA libraries were constructed by using the Illumina Truseq™ RNA sample prep Kit (Illumina, USA) following the manufacturer's instruction. The libraries were sequenced using Illumina HiSeq4000 instrument by Shanghai Majorbio Bio-pharm Technology Corporation, China. All clean Illumina sequencing data have been uploaded to SRA database.
3. Results and discussion 2.3. De novo assembly and gene annotation 3.1. Transcriptome sequencing and assembly After data filtering, poor quality reads such as adaptors, low quality sequences (< Q20), and the reads less than 20 bp in length were removed. The remaining reads were de novo assembled by software Trinity (http://trinityrnaseq.sourceforge.net/) [24]. All assembled sequences were searched against the Pfam database (http://pfam.sanger.ac.uk/) using HMMER3 program (http://hmmer.
Six cDNA libraries were generated, and a total of 175.40 and 177.52million clean reads for the wild and mutant ZI larvae were obtained, respectively (Table 1, SRA: PRJNA595668). The sequencing quality was high, with Q20 ratio > 96% and Q30 ratio > 90% of each sample. And the average GC content was 51.96% (Table 1). 516
Fish and Shellfish Immunology 98 (2020) 515–521
Y. Zhao, et al.
Table 1 Summary of transcriptome sequencing and assembly of E. sinensis. Parameters
Reads Total raw reads (Mb) Total Clean reads (Mb) Total Clean bases (Gb) Clean reads Q20 (%) Clean reads Q30 (%) GC percentage (%)
Assembly statistics Total Number Average length (bp) GC percentage (%) N50 of transcripts (bp)
Wild groups
Mutational groups
ES1
ES2
ES3
ESR1
ESR2
ESR3
64.39 61.63 8.80 96.67 90.50 52.82
55.60 53.47 7.67 96.70 90.54 53.12
62.42 60.26 8.67 96.85 90.89 52.44
59.06 56.51 8.08 96.71 90.61 51.10
59.87 57.59 8.24 96.85 90.98 51.04
65.55 63.42 9.11 97.00 91.30 51.21
transcripts 82360 1307 49.45 2249
unigenes 63387 1065 48.85 1758
After de novo assembly using Trinity, 82360 transcripts and 63387 unigenes were assembled with an average length of 1307bp and 1065bp, respectively. The GC content of these assembled transcripts was 49.45%, and N50 length was 2249 bp (Table 1). These values were higher than previous study of E. sinensis transcriptome-analysis results [30]. The most of the transcripts were less than 800 bp (51.18%), while 19.13% transcripts were over 2000 bp. These outcomes indicated that the sequences data of present study were high-quality and reliable due to its high quality, and more information on relevant genes were expected to be collected.
Fig. 2. Venn diagram of unigene numbers annotated in five databases.
single-organism process (3750, 40.6%) and regulation of biological process (1549, 16.8%) subcategories in biological process category. Besides, 66 unigenes were enriched in immune system process (0.7%) (Fig. S1). 6538 unigenes were categorized into 25 categories in KOG database (Fig. S2, Table S2). Among them, the largest number of unigenes were assigned into the categories of translation, ribosomal structure and biogenesis (915, 14.0%), followed by signal transduction mechanisms (774, 11.8%) and posttranslational modification, protein turnover and chaperones (611, 9.3%), successively. On the contrary, the smallest numbers of unigenes were assigned into extracellular structures (20, 0.3%), nuclear structure (10, 0.2%) and cell motility (5, 0.1%). In addition, 24 unigenes (0.4%) were assigned to defense mechanisms. After KEGG pathway analysis, 17132 unigenes were classified into 372 known pathways and the category with the most unigenes is Ribosome (1263, 7.3%), followed by Endocytosis (526, 3.1%), Huntington's disease (495, 2.9%), Protein processing in endoplasmic reticulum (461, 2.7%), and RNA transport (450, 2.6%). The 372 pathways were grouped into five biochemical pathways, including metabolism (A), genetic information processing (B), environmental information processing (C), cellular processes (D) and organismal systems (E) (Fig. 4). “Global and overview maps” (2671, 15.6%), was the most abundant subcategories within the category of metabolism. As for genetic information processing, environmental information processing and cellular processes, the most common subcategory was “Translation” (2382, 13.9%), “Signal transduction” (1942, 11.3%) and “Cell growth and death” (1026, 6.0%), respectively. “Endocrine system” (1221, 7.1%) and “immune system” (983, 5.7%) were the two most abundant subcategories in organismal systems (Table S3).
3.2. Functional annotation and classification Total of 63387 unigenes were annotated with mean length of 1065 bp, and N50 length of 1758 bp were further assembled. Similar with the transcripts, most unigenes were less than 800 bp (58.48%), and many transcripts were more than 2000 bp in length (13.66%). All of the unigenes were annotated in at least one database, and 3373 genes (5.32%) of them were annotated in all five databases. The annotation success rate of unigenes was 22015 in NR (34.73%), 17132 in KEEG (27.03%), 20919 in SwissProt (33.00%), 17316 in Pfam (27.32%), and 6538 in String (10.31%), respectively (Table 2, Fig. 2). The species distribution revealed that E. sinensis has the highest number of hits to the Zootermopsis nevadensis (1086, 8.2%), followed by Daphnia pulex (1290, 5.8%), Tribolium castaneum (745, 3.4%), Pediculus humanus (500, 2.3%), Nasonia vitripennis (371, 1.7%) and other species (Fig. 3). 3.3. Functional classification, enrichment, and pathway analysis 9243 of the 12236 predicted transcripts were successfully assigned to GO annotation within the three primary GO categories (molecular functional, cellular component, and biological process) and 59 sub-categories based on the NR annotation. Most unigenes were assigned to metabolic process (5273, 57.0%), cellular process (5133, 55.5%),
3.4. Identification and analysis of DEGs
Table 2 Statistics of unigene functional annotation. Annotation in Database
Unigene No.
Percentage (%)
NR KEGG Swissprot Pfam String All annotated transcripts
22015 17132 20919 17316 6538 3373
34.73 27.03 33.00 27.32 10.31 5.32
After the filtration by the fold changes, 7153 DEGs were identified between the wild and mutant ZI larvae, with 5194 of them were upregulated and 1959 were down-regulated. Statistical significance as well as magnitude of change was visualized by volcano plot (Fig. S3). The degree in change of up-regulated DEGs was apparently higher than that of down-regulated DEGs. GO function and pathway enrichment analysis divided all DEGs into three categories, biological process, cellular component and molecular 517
Fish and Shellfish Immunology 98 (2020) 515–521
Y. Zhao, et al.
Fig. 3. Comparison of E. sinensis transcriptomic sequences with the known sequences in database.
membrane-bounded organelle”, and “non-membrane-bounded organelle” (p < 0.001). KEGG pathway analysis assigned a total of 3539 DEGs into 248 known pathways, which belonged to five core categories, cellular process, metabolism, environmental information processing, genetic information processing and organismal systems. “Transport and catabolism” (223, 45.9%) in the category of cellular process, “Carbohydrate metabolism” (240, 22.4%) in the category of metabolism, “Signal transduction” (394, 87.6%) in the category of environmental
function. The highest percentage of DEGs involved in GO terms was “metabolic process” of the biological processes (613, 8.6%), “cellular component” (607, 8.5%), in the category of cellular component, and “molecular function” (930, 13.0%) in the category of molecular functions. Finally, the top ten significantly changed categories were as followed, “ribonucleoprotein complex biogenesis”, “structural constituent of ribosome”, “cellular macromolecule biosynthetic process”, “ribosome biogenesis”, “structural molecule activity”, “ribosome”, “translation”, “macromolecule biosynthetic process”, “intracellular non-
Fig. 4. KEGG functional classification of assembled unigenes. 518
Fish and Shellfish Immunology 98 (2020) 515–521
Y. Zhao, et al.
Table 3 Immune-related pathways in KEGG and Go Terms. Database KEGG ko04612 ko04640 ko04670 ko04622 ko04611 ko04621 ko04610 ko04666 ko04660 ko04062 ko04620 ko04662 ko04650 GO GO:0002697 GO:0002698 GO:0002683
Pathway
Up No.
Down No.
Total No.
p value
Antigen processing and presentation Hematopoietic cell lineage Leukocyte transendothelial migration RIG-I-like receptor signaling pathway Platelet activation NOD-like receptor signaling pathway Complement and coagulation cascades Fc gamma R-mediated phagocytosis T cell receptor signaling pathway Chemokine signaling pathway Toll-like receptor signaling pathway B cell receptor signaling pathway Natural killer cell mediated cytotoxicity
33 2 10 1 13 6 1 14 3 11 0 2 1
6 4 8 3 8 2 1 14 2 11 3 0 1
39 6 18 4 21 8 2 28 5 22 2 2 2
0.00501 0.28483 0.46669 0.51429 0.52441 0.63878 0.70776 0.92131 0.96407 0.99027 0.99087 0.99449 0.99821
regulation of immune effector process negative regulation of immune effector process negative regulation of immune system process
2 1 1
1 1 1
3 2 2
0.00526 0.0127 0.0353
information processing, the “Translation” (386, 66.3%) subcategory of “genetic information processing” category, and “Cardiac muscle contraction” (53) of organismal systems category was the most common subcategories, respectively. Moreover, the most abundant pathways included “Ribosome” (293, 23.2%), “Oxidative phosphorylation” (89, 22.6%), “Lysosome” (76, 21.3%) and “Carbon metabolism” (67, 15.0%) and “Protein processing in endoplasmic reticulum” (66, 14.3%). The top five significantly distinct categories were “Ribosome”, “Oxidative phosphorylation”, “Lysosome”, “Cardiac muscle contraction” and “Proteasome” (p < 0.001).
ATPase subunit alpha (ATP), glutamate receptor 1 (AMPAR), glutamate receptor (NMDAR), cAMP-dependent protein kinase catalytic subunit (PKA), Ras-related protein Rap-1 (Rap1), Ras-like GTP-binding protein RHO (Rho), serine/threonine-protein phosphatase PP1 catalytic subunit (PP1), Aconitate hydratase (ACO), p21-activated kinase 1 (PAK1), phosphatidylinositol phospholipase C epsilon (PLCε) and Ras homolog gene family member A (RhoA) in the cAMP signaling pathway, as well as low density lipoprotein receptor-related protein 5/6 (LRP5/6), casein kinase II subunit alpha (CK2), glycogen synthase kinase 3 beta (GSK3β), PKA, RhoA, serine/threonine-protein phosphatase 2B catalytic subunit (CaN) and C-terminal binding protein (CtBP) in the Wnt signaling pathway. Most of those DEGs exhibited significantly up-regulated expression. It is well known that melanin has various functions in different organisms, such as protecting cell from UV light damage for human, improving avian feathers’ resistance to abrasion, etc. [42]. Beyond that, for crustacean, melanin plays a significant role in immune response by involving in wound melanization, encapsulation of foreign particles, clot formation, and so on [43,44]. Although the underlying mechanisms are still unclear, those DEGs may make a significant contribution to the shell color variation of E. sinensis and the immune deficiencies of brownish-orange E. sinensis. Moreover, GO analysis identified 7 DEGs enriched in the “Immune system process” including “regulation of immune effector process”, “negative regulation of immune effector process” and “negative regulation of immune system process” (p < 0.05) (Table 3). Specifically, one of lipocalins, Maltose-binding periplasmic protein (MBP), was enriched in all the three processes. As transport proteins, lipocalin family takes part in many biological functions, including retinol transport, pheromone transport, prostaglandin synthesis, invertebrate camouflage, olfaction, immune response, etc. [45]. In the present study, 34 DEGs were related to the lipocalins, and 13 of them were significantly different, including the Apolipoprotein D-like, crustacyanin-like lipocalin, putative cellular retinoic acid/retinol binding protein (CRABP-II), MBP, Fatty acid-binding protein and Apolipoprotein D (ApoD). Among them, the Apolipoprotein D-like, crustacyaninlike lipocalin, and MBP were all down-regulated whereas CRABP-II, Fatty acid-binding protein and ApoD were up-regulated in mutant E. sinesis. This may illustrate that up/down-regulation of lipocalin family proteins were related to body color of E. sinesis.
3.5. Analysis of DEGs of immune-related genes Finally, 13 KEGG pathways related to immune system were detected among 248 pathways mapped to the KEGG database. The numbers of DEGs in the 13 pathways were listed in Table 3. Our KEGG analysis showed that the up/down-regulated DEGs were present in 13 “immune system” pathways, and enrichment in most pathways was not significant except pathway “Antigen processing and presentation” (p < 0.05) (Table 3). 33 up-regulated and 6 down-regulated DEGs were enriched on ten unigenes, including HSP70 (heat shock protein 70), HSP90 (heat shock protein 90), CTSB/L/S (cathepsin B/L/S) and so on (Fig. S4). HSP can protect crustacean from various internal or external stresses, and it is evidently involved in innate immune response [31–34]. CTSB was reported to participate in the regulation of developmental processes in some arthropods [35]. In addition, as an important participant in the MHC class I pathway, most of the proteasome subunits are significantly up-regulated in mutant type (Fig. S5). The proper expression of proteasome plays an important role in cellular homeostasis, and proteasome dysregulation have been associated with many human diseases, such as cancer and Parkinson's disease [36–39]. Similarly, abnormal proteasome expression was also found in mud crab and silkworm in their response to viral or fungal infection [40,41]. In the present study, combined with the overexpression of heat shock protein and proteasome, it may infer that the mutant individuals are experiencing some kind of stress (probably caused by genetic defect). In this study, with the exception of pathway “Antigen processing and presentation”, some significantly differential DEGs have been located in melanogenesis related pathways, such as cyclic adenosine monophosphate/protein kinase A (cAMP/PKA) signaling pathway, Wnt signaling pathway, mitogen-activated protein kinases (MAPK) signaling pathway, and so on. For example, we found DEGs including adrenergic receptor beta-1(GPCR), adenylate cyclase 1 (AC), cyclic nucleotide gated channel alpha 1 (CNGC), Calmodulin (CaM), Plasma membrane calcium-transporting ATPase (PMCA), sodium/potassium-transporting
3.6. Validation of ten DEGs by qRT-PCR The expression of ten chosen DEGs were verified on three individual biological replicates through qRT-PCR. Among ten chosen genes, four of them were down-regulated in mutant type (Astacin, Lipocalin, MBP 519
Fish and Shellfish Immunology 98 (2020) 515–521
Y. Zhao, et al.
Appendix A. Supplementary data Supplementary data to this article can be found online at https:// doi.org/10.1016/j.fsi.2020.01.051. References [1] N.G. Ertl, A. Elizur, P. Brooks, A.V. Kuballa, T.A. Anderson, W.R. Knibb, Molecular characterisation of colour formation in the prawn Fenneropenaeus merguiensis, PloS One 8 (2013) e56920, , https://doi.org/10.1371/journal.pone.0056920. [2] Y. Murakami, K. Wada, Inter-populational variations in body color related to growth stage and sex in Gaetice depressus (De Haan, 1835) (Decapoda, Brachyura, Varunidae), Crustaceana 88 (2015) 113–126, https://doi.org/10.1163/1568540300003391. [3] A.S. Freire, M.A.A. Pinheiro, H. Karam-Silva, M.M. Teschima, Biology of Grapsus grapsus (Linnaeus, 1758) (Brachyura, Grapsidae) in the Saint Peter and Saint Paul Archipelago, Equatorial Atlantic ocean, Helgol. Mar. Res. 65 (2011) 263–273, https://doi.org/10.1007/s10152-010-0220-5. [4] M.Z. Darnell, Ecological physiology of the circadian pigmentation rhythm in the fiddler crab Uca panacea, J. Exp. Mar. Biol. Ecol. 426–427 (2012) 39–47, https:// doi.org/10.1016/j.jembe.2012.05.014. [5] S.M. Kronstadt, M.Z. Darnell, P. Munguia, Background and temperature effects on Uca panacea color change, Mar. Biol. 160 (2013) 1373–1381, https://doi.org/10. 1007/s00227-013-2189-5. [6] R.K. Tume, A.L. Sikes, S. Tabrett, D.M. Smith, Effect of background colour on the distribution of astaxanthin in black tiger prawn (Penaeus monodon): effective method for improvement of cooked colour, Aquaculture 296 (2009) 129–135, https://doi.org/10.1016/j.aquaculture.2009.08.006. [7] M. Ikeda, T. Suzuki, Y. Fujio, Genetic differentiation among populations of Japanese freshwater crab, Geotheiphusa dehaani (White), with reference to the body color variation, Benthos Res. 53 (1998) 47–52, https://doi.org/10.5179/benthos1996. 53.1_47. [8] N.H. Nguyen, J. Quinn, D. Powell, A. Elizur, N.P. Thoa, J. Nocillado, R. Lamont, C. Remilton, W. Knibb, Heritability for body colour and its genetic association with morphometric traits in Banana shrimp (Fenneropenaeus merguiensis), BMC Genet. 15 (2014) 1–12, https://doi.org/10.1186/s12863-014-0132-5. [9] J.J. Nordlund, R.E. Boissy, V.J. Hearing, R.A. King, W.S. Oetting, J.P. Ortonne, The Pigmentary System: Physiology and Pathophysiology, second ed., Wiley–Blackwell, 2007 1–1229. [10] W. Zhu, L. Wang, Z. Dong, X. Chen, F. Song, N. Liu, H. Yang, J. Fu, Comparative transcriptome analysis identifies candidate genes related to skin color differentiation in red tilapia, Sci. Rep. 6 (2016) 1–12, https://doi.org/10.1038/srep31347. [11] D. Ma, H. Yang, L. Sun, D. Xu, Comparative analysis of transcriptomes from albino and control sea cucumbers, Apostichopus japonicus, Acta Oceanol. Sin. 33 (2014) 55–61, https://doi.org/10.1007/s13131-014-0464-z. [12] L. Xing, L. Sun, S. Liu, X. Li, L. Zhang, H. Yang, Transcriptome analysis provides insights into the mechanism of albinism during different pigmentation stages of the albino sea cucumber Apostichopus japonicus, Aquaculture 486 (2018) 148–160, https://doi.org/10.1016/j.aquaculture.2017.12.016. [13] D. Feng, Q. Li, H. Yu, X. Zhao, L. Kong, Comparative transcriptome analysis of the pacific oyster Crassostrea gigas characterized by shell colors: identification of genetic bases potentially involved in pigmentation, PloS One 10 (2015) 1–17, https://doi. org/10.1371/journal.pone.0145257. [14] A.T. Palma, R.S. Steneck, Does variable coloration in juvenile marine crabs reduce risk of visual predation ? Ecology 82 (2001) 2961–2967, https://doi.org/10.2307/ 2679974. [15] M.F. Tlusty, A. Metzler, S. Huckabone, S. Suanda, S. Guerrier, Morphological colour change in the American lobster (Homarus americanus) in response to background colour and UV light, N. Z. J. Mar. Freshw. Res. 43 (2009) 247–255, https://doi.org/ 10.1080/00288330909509998. [16] N.M. Wade, S. Cheers, N. Bourne, S. Irvin, D. Blyth, B.D. Glencross, Dietary astaxanthin levels affect colour, growth, carotenoid digestibility and the accumulation of specific carotenoid esters in the Giant Tiger Shrimp, Penaeus monodon, Aquacult. Res. 48 (2017) 395–406, https://doi.org/10.1111/are.12888. [17] M. Wang, X. Zhu, J. Yang, Z. Dai, K. Mahmood, F. Yang, W. Yang, Prawn lipocalin: characteristics and expressional pattern in subepidermal adipose tissue during reproductive molting cycle, Comp. Biochem. Physiol. B Biochem. Mol. Biol. 147 (2007) 222–229, https://doi.org/10.1016/j.cbpb.2007.01.003. [18] N.M. Wade, A. Tollenaere, M.R. Hall, B.M. Degnan, Evolution of a novel carotenoidbinding protein responsible for crustacean shell color, Mol. Biol. Evol. 26 (2009) 1851–1864, https://doi.org/10.1093/molbev/msp092. [19] F. Yang, M. Wang, Y. Ma, W. Ma, W. Yang, Prawn lipocalin: characterization of a color shift induced by gene knockdown and ligand binding assay, J. Exp. Zool. Part A Ecol. Genet. Physiol. 315 A (2011) 562–571, https://doi.org/10.1002/jez.706. [20] D. Chen, M. Zhang, S. Shrestha, Compositional characteristics and nutritional quality of Chinese mitten crab (Eriocheir sinensis), Food Chem. 103 (2007) 1343–1349, https://doi.org/10.1016/j.foodchem.2006.10.047. [21] Y. Si, X. Li, Y. Jiang, Y. Wang, X. Liu, Y. Zheng, J. Suo, N. Sun, Culture pattern of juvenile Chinese mitten handed crab Eriocheir sinensis with red color shell, Journal of Dalian Ocean University 32 (2017) 139–144, https://doi.org/10.16535/j.cnki. dlhyxb.2017.02.003. [22] M.L. Metzker, Sequencing technologies the next generation, Nat. Rev. Genet. 11 (2010) 31–46, https://doi.org/10.1038/nrg2626. [23] M. Garber, M.G. Grabherr, M. Guttman, C. Trapnell, Computational methods for
Fig. 5. qRT-PCR verification of differentially expressed genes between the wild and mutant E. sinensis. *Significant level (p < 0.05). ** (p < 0.01),*** (p < 0.001).
and CHH1), whereas other six were up-regulated (CaM 1, CTSL, HSP70, LSH-R, MEP1A, PKA C) (Fig. 5). The results of qRT-PCR were consistent with expression pattern of RNA-Seq analyses, which support the accuracy and reliability of DEGs given by RNA-Seq. 4. Conclusion In conclusion, we constructed and sequenced the color-specific transcriptomes ZI larvae of E. sinensis. Illumina paired-end RNA sequencing technology was used to sequence. In total, More than 350 million clean reads were obtained, which then assembled into 82360 transcripts and 63387 unigenes. The unigenes were annotated by six databases, NR, Pfam, KOG, Swiss-prot, KEGG, and GO, and we obtained comprehensive gene function information. We identified 7153 significantly DEGs including 5194 significantly up-regulated and 1959 down-regulated DEGs. KEGG and GO analysis indicated those DEGs in pathway “Antigen processing and presentation”, melanogenesis related pathway, as well as those of lipocalin family protein may associated with shell color alteration and immunodeficiency of mutant E. sinensis. Overall, this study provides valuable information on the molecular mechanisms underlying interactions between color determination and the immune system in E. sinensis, but also facilitates further investigations of functional genomics for this species and other closely related species. Author contributions Y Zhao and X Li conceived and designed the experiments. Z Han, Y Zhang, T Dong, H Wei, J Dong and Y Li performed the experiments. Y Zhao and X Zhu analyzed the data and wrote the paper. X Li contributed the reagents and materials. Acknowledgments This work was supported by grants from the Cultivation Plan for Youth Agricultural Science and Technology Innovative Talents of Liaoning Province (No. 2015044), Talent Introduction Program of Shenyang Agricultural University (No. 880416005), Guide Project for Key Research and Development Program of Liaoning province (No. 201802120), and Science and Technology Plan Project of Liaoning Province (No. 2019JH2/10200006). 520
Fish and Shellfish Immunology 98 (2020) 515–521
Y. Zhao, et al.
[24]
[25]
[26]
[27]
[28]
[29]
[30]
[31]
[32]
[33]
237–246, https://doi.org/10.1016/j.fsi.2016.11.049. [34] X. Chu, T. Yang, Y. Liu, L. Hong, T. Jiao, X. Meng, D. Zhang, J. Wang, B. Tang, C. Zhou, Q. Liu, W. Zhang, W. He, Transcriptome analysis of differential expressed genes in hepatopancreas of Procambarus clarkii challenged with peptidoglycan, Fish Shellfish Immunol. 86 (2019) 311–318, https://doi.org/10.1016/j.fsi.2018.11.048. [35] H. Wang, C. Fu, Q. Zeng, F. Li, H. Wang, J. Sun, Comparative transcriptome analysis reveals related regulatory mechanisms of androgenic gland in Eriocheir sinensis, BioMed Res. Int. 1 (2017) 1–12, https://doi.org/10.1155/2017/4956216. [36] A. Kumatori, K. Tanaka, N. Inamura, S. Sone, T. Ogura, T. Matsumoto, T. Tachikawa, S. Shin, A. Ichihara, Abnormally high expression of proteasomes in human leukemic cells, Proc. Natl. Acad. Sci. U. S. A 87 (1990) 7071–7075, https:// doi.org/10.1073/pnas.87.18.7071. [37] L. Chen, K. Madura, Increased proteasome activity, ubiquitin-conjugating enzymes, and eEF1A translation factor detected in breast cancer tissue, Canc. Res. 65 (2005) 5599–5606, https://doi.org/10.1158/0008-5472.CAN-05-0201. [38] Q. Zheng, T. Huang, L. Zhang, Y. Zhou, H. Luo, H. Xu, X. Wang, Dysregulation of ubiquitin-proteasome system in neurodegenerative diseases, Front. Aging Neurosci. 8 (2016) 1–10, https://doi.org/10.3389/fnagi.2016.00303. [39] R. Motosugi, S. Murata, Dynamic regulation of proteasome expression, Front. Mol. Biosci. 6 (2019) 4–11, https://doi.org/10.3389/fmolb.2019.00030. [40] C. Hou, G. Qin, T. Liu, T. Geng, K. Gao, Z. Pan, H. Qian, X. Guo, Transcriptome analysis of silkworm, Bombyx mori, during early response to Beauveria bassiana challenges, PloS One 9 (2014) 1–13, https://doi.org/10.1371/journal.pone. 0091189. [41] S. Liu, G. Chen, H. Xu, W. Zou, W. Yan, Q. Wang, H. Deng, H. Zhang, G. Yu, J. He, S. Weng, Transcriptome analysis of mud crab (Scylla paramamosain) gills in response to Mud crab reovirus (MCRV), Fish Shellfish Immunol. 60 (2017) 545–553, https://doi.org/10.1016/j.fsi.2016.07.033. [42] F. Solano, Melanins: skin pigments and much more—types, structural models, biological functions, and formation routes, New J. Sci. 2014 (2014) 1–28, https:// doi.org/10.1155/2014/498276. [43] H. Bilandžija, M. Laslo, M.L. Porter, D.W. Fong, Melanization in response to wounding is ancestral in arthropods and conserved in albino cave species, Sci. Rep. 7 (2017) 1–11, https://doi.org/10.1038/s41598-017-17471-2. [44] L. Vazquez, J. Alpuche, G. Maldonado, C. Agundis, A. Pereyra-Morales, E. Zenteno, Immunity mechanisms in crustaceans, Innate Immun. 15 (2009) 179–188, https:// doi.org/10.1177/1753425909102876. [45] D.R. Flower, The lipocalin protein family: structure and function, Biochem. J. 318 (1996) 1–14, https://doi.org/10.1042/bj3180001.
transcriptome annotation and quantification using RNA-seq, Nat. Methods 8 (2011) 469–477, https://doi.org/10.1038/nmeth.1613. M.G. Grabherr, B.J. Haas, M. Yassour, J.Z. Levin, D.A. Thompson, I. Amit, X. Adiconis, L. Fan, R. Raychowdhury, Q. Zeng, Z. Chen, E. Mauceli, N. Hacohen, A. Gnirke, N. Rhind, F. Di Palma, B.W. Birren, C. Nusbaum, K. Lindblad-Toh, N. Friedman, A. Regev, Full-length transcriptome assembly from RNA-Seq data without a reference genome, Nat. Biotechnol. 29 (2011) 644–652, https://doi.org/ 10.1038/nbt.1883. A. Conesa, S. Götz, J.M. García-Gómez, J. Terol, M. Talón, M. Robles, Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research, Bioinformatics 21 (2005) 3674–3676, https://doi.org/10.1093/ bioinformatics/bti610. B. Li, C.N. Dewey, RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome, BMC Bioinf. 12 (2011) 323–338, https://doi. org/10.1093/10.1186/1471-2105-12-323. M.D. Robinson, D.J. McCarthy, G.K. Smyth, edgeR: a Bioconductor package for differential expression analysis of digital gene expression data, Bioinformatics 26 (2009) 139–140, https://doi.org/10.1093/bioinformatics/btp616. D.V. Klopfenstein, L. Zhang, B.S. Pedersen, F. Ramírez, A.W. Vesztrocy, A. Naldi, C.J. Mungall, J.M. Yunes, O. Botvinnik, M. Weigel, W. Dampier, C. Dessimoz, P. Flick, H. Tang, GOATOOLS: a Python library for Gene Ontology analyses, Sci. Rep. 8 (2018) 1–17, https://doi.org/10.1038/s41598-018-28948-z. C. Xie, X. Mao, J. Huang, Y. Ding, J. Wu, S. Dong, L. Kong, G. Gao, C. Li, L. Wei, KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases, Nucleic Acids Res. 39 (2011) 316–322, https://doi.org/10.1093/nar/ gkr483. Y. Liu, M. Hui, Z. Cui, D. Luo, C. Song, Y. Li, L. Liu, Comparative transcriptome analysis reveals sex-biased gene expression in juvenile Chinese mitten crab Eriocheir Sinensis, PloS One 10 (2015) 1–21, https://doi.org/10.1371/journal.pone.0133068. L. Cerenius, P. Jiravanichpaisal, H. Liu, I. Söderhäll, Crustacean Immunity, Invertebrate Immunity, Springer, Boston, MA, 2010, pp. 239–259, https://doi.org/ 10.1007/978-1-4419-8059-5_13. M.K. Chaurasia, F. Nizam, G. Ravichandran, M.V. Arasu, N.A. Al-Dhabi, A. Arshad, P. Elumalai, J. Arockiaraj, Molecular importance of prawn large heat shock proteins 60, 70 and 90, Fish Shellfish Immunol. 48 (2016) 228–238, https://doi.org/10. 1016/j.fsi.2015.11.034. W. Junprung, P. Supungul, A. Tassanakajon, HSP70 and HSP90 are involved in shrimp Penaeus vannamei tolerance to AHPND-causing strain of Vibrio parahaemolyticus after non-lethal heat shock, Fish Shellfish Immunol. 60 (2017)
521