Hepatopancreas immune response during molt cycle in the mud crab, Scylla paramamosain

Molt is a critical developmental process in crustaceans. Recent studies have shown that the hepatopancreas is an important source of innate immune molecules, yet hepatopancreatic patterns of gene expression during the molt cycle which may underlie changes in immune mechanism are unknown. In this study, we performed Illumina sequencing for the hepatopancreas of the mud crab, Scylla paramamosain during molt cycle (pre-molt stage, post-molt stage, and inter-molt stage). A total of 44.55 Gb high-quality reads were obtained from the normalized cDNA of hepatopancreas. A total of 70,591 transcripts were assembled; 55,167 unigenes were identified. Transcriptomic comparison revealed 948 differentially expressed genes (DEGs) in the hepatopancreas from the three molt stages. We found that genes associated with immune response patterns changed in expression during the molt cycle. Antimicrobial peptide genes, inflammatory response genes, Toll signaling pathway factors, the phenoloxidase system, antioxidant enzymes, metal-binding proteins and other immune related genes are significantly up-regulated at the post-molt stage and inter-molt stage compared with the pre-molt stage, respectively. These genes are either not expressed or are expressed at low levels at the pre-molt stage. To our knowledge, this is the first systematic transcriptome analysis of genes capable of mobilizing a hepatopancreas immune response during the molt cycle in crustaceans, and this study will contribute to a better understanding of the hepatopancreas immune system and mud crab prophylactic immune mechanisms at the post-molt stage.

In crustaceans, molt is an essential characteristic of the organism's development; in order to grow, the animal must shed the old exoskeleton and construct a new one [1][2][3][4] . The crustacean molt cycle affects the status of a number of physiological processes, including the interaction with environmental stressors and disease agents 5,6 . In crustaceans, several researchers have described studies of molt susceptibility to pathogenic infections. Corteel et al. 7 reported that the Pacific whiteleg shrimp, Litopenaeus vannamei at post-molt stage is more susceptible to white spot syndrome virus (WSSV) infection than at pre-molt stage. Tumburu et al. 8 have also shown that pesticide-virus interactions lead to increasingly higher susceptibility to acute toxicity at post-molt stage in L. vannamei. These studies indicate that crustaceans display susceptibility to pathogenic infections that vary during the molt cycle.
The hepatopancreas is an important organ involved in the process of crustacean molting, and plays a vital role in energy storage and breakdown, nutrient accumulation, and carbohydrate and lipid metabolism 9 . Growth and developmental metabolism related genes that are involved in the formation of the molt have been isolated in crustaceans. However, despite extensive research, the changes in hepatopancreas immune response during the molt cycle in crustaceans still remain poorly understood. Recent findings show that the hepatopancreas is an integrated organ of immunity and metabolism, which is an important source of innate immune molecules in crustaceans [10][11][12][13] , while nutrition and metabolic regulation may have an impact on the efficiency of the immune response in the hepatopancreas 14 . At post-molt stage, in light of the soft body, lack of exoskeleton protection, and susceptibility to various pathogens, crustaceans may mobilize the body's immune defense 4 .
Mud crabs (Scylla species) are found throughout tropical and warm temperate zones in the Indo-Pacific, and S. paramamosain is a commercially important marine fishery species in the Southeast coast of China 15 . Transcriptomics is an emerging discipline, which can now be used to examine patterns of gene transcription relative to crustacean immunity, growth and molting 12,14,16 . In this study, we used RNA-Seq to investigate the Functional classification of DEGs of hepatopancreas. Assembled non-redundant unigenes were subjected to Gene Ontology (GO) 17 , Clusters of Orthologous Groups (COG) 18 and Kyoto Encyclopedia of Genes and Genomes (KEGG) 19 databases for blast searching 20 . GO annotation output showed that the unigenes could be assigned to three parts: biological processes, cellular components, and molecular functions. GO function classification was performed with DEGs. In the hepatopancreas, compared post-molt stage and inter-molt stage, the three major functional categories (biological process, cellular component, and molecular function) were not found in the GO terms associated with significant DEGs. When inter-molt stage and pre-molt stage were compared, the DEGs were assigned to 36 terms, of which dominant subcategories included single-organism process (64 genes) and metabolic process (84 genes) in the category of biological process, organelle (23 genes) and cell part (23 genes) in cellular component, and catalytic activity (149 genes) in molecular function, respectively (Fig. 1a). In the comparison of pre-molt stage and post-molt stage, metabolic process (33 genes) was the top abundant terms in biological process. While, extracellular region (6 genes) and catalytic activity (49 genes) were the top terms in the cellular component and molecular function, respectively (Fig. 1b) (Supplementary data 1).  www.nature.com/scientificreports/ To classify orthologous gene products, COG function classification was performed with significant DEGs. In the hepatopancreas, DEGs at post-molt stage and inter-molt stage were not found in the COG annotations. In the comparison of inter-molt stage and pre-molt stage, the cluster of "general function prediction only" (36,16.51%) represented the largest group, followed by "amino acid transport and metabolism" (34,15.6%), "carbohydrate transport and metabolism" (31, 14.22%), "Secondary metabolism biosynthesis, transport and catabolism" (26,11.9%), "lipid transport and metabolism" (25,11.47%) and "inorganic ion transport and metabolism" (15, 6.88%), whereas "coenzyme transport and metabolism" (1, 0.46%) was the group with the lowest number of identifiable consensus sequences (Fig. 2a). In the comparison of pre-molt stage and post-molt stage, the cluster of "carbohydrate transport and metabolism" (15,23.81%) represented the largest group, followed by "amino acid transport and metabolism" (10, 15.87%), "lipid transport and metabolism" (9, 14.29%), "general function prediction only" (9, 14.29%), "Secondary metabolism biosynthesis, transport and catabolism" (8, 12. 7%) and "energy production and conversion" (4, 6.38%), whereas "coenzyme transport and metabolism" (1, 1.59%) was the group with the lowest number of identifiable consensus sequences (Fig. 2b) (Supplementary data 2).
KEGG enrichment results were not obtained from 6 DEGs in the post-molt stage and inter-molt stage comparisons. KEGG enrichment was carried out on 789 DEGs obtained from an inter-molt stage and pre-molt stage comparison, among which 268 DEGs were classified into cellular processes, environmental information processing, metabolism (Fig. 3a). In cellular processes, fourteen genes were assigned to lysosome (ko: 04142) and four genes were assigned to phagosome (ko: 04145); In environmental information processing, one gene related to signal transduction was assigned to FoxO signaling pathway (ko: 04068); In metabolism, carbon metabolism (ko: 01200) and glycine, serine and threonine metabolism (ko: 00260) were the top two pathway (21genes). KEGG enrichment was carried out on 153 DEGs obtained from a pre-molt stage and post-molt stage comparison, among which 67 DEGs were classified into cellular processes, metabolism and organismal systems (Fig. 3b). In cellular processes, five genes were assigned to lysosome (ko: 04142) and one gene was assigned to phagosome (ko: 04145); In metabolism, amino sugar and nucleotide sugar metabolism (ko: 00520) and starch and sucrose metabolism (ko: 00500) were the top two pathways (14 genes); In organismal systems, one gene related to carbohydrate digestion and absorption was identified (ko: 04973) (Supplementary data 3).
DEGs in the hepatopancreas transcriptome related to immune system. To date, there have been few studies on the immune system in the hepatopancreas of the mud crab. DEGs identified across the different molt stages (post-molt stage, inter-molt stage, and pre-molt stage) were compared pair-wise, and genes related to immune response were screened. Expression abundance of the genes showed a significant difference. Among them, we focused on key genes involved in Toll, immune deficiency (IMD) signaling pathways, immune-related and cytokines (Table 3). We found a large number of highly significant up-regulated immune response genes during the post-molt stage and the inter-molt stage, however, these genes are significantly down regulated and not expressed at the pre-molt stage ( transcript validation by qRt-pcR. Twelve immune response genes that identified as differentially expressed genes (DEGs) at post-molt and inter-molt relative to pre-molt were examined by qRT-PCR, and the expression levels of these immune genes were significantly associated with the RNA-seq results (Fig. 5). Our results confirmed the reliability of RNA-seq and accuracy of the Trinity assembly.

Discussion
In this study, we assembled the hepatopancreas transcriptomes at post-molt, inter-molt and pre-molt stages, and obtained an annotated transcriptome over the molt cycle of the mud crab, S. paramamosain. We found classifications of DEGs from GO, COG and KEGG databases, which were involved in energy metabolism, growth development, osmotic pressure regulation and immune-related processes, this phenomenon demonstrated the major physiological changes happen during molt cycle in S. paramamosain. Compared with pre-molt stage, a large number of genes were up-regulated at the post-molt stage and the inter-molt stage, which further showed that the hepatopancreas is an important organ involved in molt cycle in S. paramamosain. It is interesting to note that some immune-related genes (antimicrobial peptide genes, inflammatory response genes, Toll signaling pathway factors, phenoloxidase system, antioxidant enzymes, metalloproteins, etc.) were significantly up-regulated at the post-molt stage and the inter-molt stage compared with the pre-molt stage, respectively. This might indicate that the mud crabs are sensitive to environmental and physiological changes, and therefore adopt prophylactic immune mechanisms. Similar patterns exist in the life cycle of the swimming crab, Portunus trituberculatus, the chitinase gene PtCht-1 involved in immune defense is highly expressed at post-molt stage and inter-molt stage, but has a low expression at pre-molt stage 21 . It was also reported that the expression of antimicrobial peptides (AMPs) during molt cycle is associated with prophylactic innate immunity in the fruit fly Drosophila melanogaster, which acts to prevent infection 22 . Antimicrobial proteins (AMps). Three unique AMP genes were discovered in the hepatopancreas cDNA library, and they showed a significant upregulation at post-molt stage and inter-molt stage.
Crustins are a cysteine-rich cationic peptide with a characteristic four-disulphide core-containing whey acidic protein (WAP) domain 23 , which has antimicrobial activity against Gram-positive bacteria and no activity against Gram-negative bacteria 24 .
SpHyastatin was first screened as a new cationic antimicrobial peptide in the suppression subtractive hybridization (SSH) cDNA library constructed from the hemocytes of S. paramamosain 25,26 . The previous study has proved that SpHyastatin can effectively kill microbes via a permeabilizing mechanism and bind to various anionic cell wall components against different species of bacteria 27  www.nature.com/scientificreports/ Lysozymes are immune effectors widely distributed and involved in many physiological processes, such as in immune and digestive systems 28,29 , exerting cytosolic activity on peptidoglycans of bacterial cell walls to initiate cell lysis 30 . In S. paramamosain, lysozyme exhibited moderate antimicrobial activity and possessed weak isopeptidase activity and had a strong binding activity to lipopolysaccharide from Escherichia coli and lipoteichoic acid and peptidoglycan (PGN) from Staphylococcus aureus 30 . Inflammatory gene IL-16. As a pro-inflammatory cytokine, Interleukin-16 (IL-16) was first identified as a lymphocyte chemoattractant factor in 1982 31 . IL-16 regulates migration, proliferation and activation of various cells, as well as the release of pro-inflammatory cytokines [32][33][34] . A recent study revealed that IL-16 might play a vital role in initiating innate immune responses against pathogen infection in S. paramamosain 35 . prophenol/phenol oxidase cascade. In the present study, proteases and inhibitors in the phenol oxidase system were significantly up-regulated at the post-molt stage. In crustaceans, the prophenoloxidase activating system (proPO-system) is an important immune defense system. Serine proteases (SPs) constitute one of the largest families of enzymes in the animal kingdom and play important roles in immune responses 36 . In crustaceans, trypsin-like serine proteases have been isolated from the hepatopancreas of the redclaw crayfish, Cherax quadricarinatus and the Chinese shrimp, Fenneropenaeus chinensis, which are involved in the innate immune defense against pathogens 36,37 . Serine proteinase homologs were also obtained from S. paramamosain, it is suggested that Sp-SPH protein could bind to a number of bacteria and play a key role in host defense against microbe invasion 38 .
Pacifastin and alpha2-macroglobulin (α2M), ubiquitous protease inhibitors, are involved in the regulation of the prophenyloxidase (proPO) system to avoid the deleterious effects of its active components [39][40][41] . The expression of these genes is presumed to be due to the activation of the proPO serine proteinase cascade, the over-activated protease is very harmful to cells and tissues, so the expression of serine protease inhibitors is used to control excessive immune responsiveness to maintain homeostasis 39,40 . toll signaling pathway. In this study, the expression of the gene coding for the Toll-like receptor was significantly up-regulated at the post-molt stage and inter-molt stage. The Toll signaling pathway is one of the most important innate immune pathways in invertebrates and has a homologous similarity to the mammalian TLR pathway 41 . Toll/Toll-like receptors (TLRs) were the first pattern recognition receptors (PRRs) identified to play a crucial role in innate immune responses in crustaceans. Up to now, TLRs have been cloned and characterized in a limited number of crustacean species, which suggested that TLRs could regulate the transcripts of several AMPs (SpALF1-6, SpCrustin, SpHistin, SpArasin, SpGRP and SpHyastatin) [42][43][44][45][46][47] . It is possible that the Toll signaling pathway is activated to regulate the expression of downstream antimicrobial peptides at the post-molt stage, and the expression of antimicrobial peptides was up-regulated in response to different potential pathogens as mentioned above. In addition, we also screened genes related to the Toll pathway, such as Toll-6, dorsal, pellelike kinase, Myd88 and tumor necrosis factor receptor-associated factor 6 (TRAF6). According to the FPKM data, there are no significant differences at the post-molt stage and inter-molt stage compared with the pre-molt stage, respectively, so they were not listed in Table 3. These genes are the homologs of Drosophila Toll pathway genes and involved in Toll-mediated immune defense 48 . It has been reported that the expression profiles of genes involved in the Toll pathway up-regulated at different time points in crustaceans 48-50 . Antioxidant enzymes. The present study found that the expression of catalase and thioredoxin were significantly up-regulated after molt, possibly due to the large amount of reactive oxygen species produced by oxidative stress at the post-molt stage. Antioxidant enzymes play essential roles in antioxidant responses caused by pathogen invasion or metabolic process 51 . Catalase is involved in the immune response and plays a role in the protection against oxidative stress in S. paramamosain 52 . Thioredoxin, with a redox active disulfide bridge, is responsible for maintaining the balance of reactive oxygen species, which has an important impact on the immune system 53,54 . Thioredoxin has been cloned and characterized in S. paramamosain and may prove to be a potential biomarker gene for evaluation of environmental stress in marine species 55 . Metal-binding proteins. Two meta-binding proteins were significantly up-regulated at the post-molt stage in S. paramamosain. We infer that the tight regulation of ferritin and zinc protease were a primary defense mechanism to resist microbial infection. Metal-binding proteins are involved in many physiological processes, ranging from immune responses to cellular signaling pathways 56 . Ferritin is a metal-chelator protein involved Figure 1. Distribution of Gene Ontology (GO) functional categories for S. paramamosain hepatopancreas transcriptome. (a) Functional classification of DEG unigenes compared between inter-molt stage and pre-molt stage in the GO database; (b) Functional classification of DEG unigenes compared between pre-molt stage and post-molt stage in the GO database; GO functional categories is mined here from https ://www.geneo ntolo gy.org. Note: the x-axis indicates the secondary classification of the GO database; the y-axis on the left side indicates the ratio of the number of genes annotated to this GO classification to all genes (the denominator of ' All Unigene' is all unigenes, and the denominator of 'DEG Unigene' is all DEG unigenes). The y-axis on the right side indicates the number of genes annotated to this GO entry. The number above represents the number of DEG unigenes, and the number below represents the number of all unigenes. www.nature.com/scientificreports/ in iron storage, and plays a crucial role in iron metabolism. As some bacteria need to isolate iron, ferritin can protect organisms from bacterial infection by regulating free iron availability in the environment 57 . Zinc protease belongs to the metal endopeptidases and is involved in a variety of physiological activities, including connective tissue remodeling and resection of nascent protein signal peptides 56 . A full-length ferritin cDNA clone has been isolated from the river prawn, Macrobrachium nipponense 58 , however, zinc protease in S. paramamosain has not yet been reported. In L. vannamei and the kuruma shrimp, Penaeus japonicas, mRNA of www.nature.com/scientificreports/ zinc proteinase and ferritin are more expressed in WSSV resistant individuals compared with susceptible ones 59,60 . These data support the role that these two metal binding proteins may guard against infection events that may accompany cuticle loss and regrowth during the molting process.
other immune factors. In the present study, macrophage migration inhibitory factor (MIF), cathepsin L (CTSL), c-type lectins (CTLs) and hemocyanin were significantly up-regulated at the post-molt stage. MIF is a regulator of the innate immune system, and plays a key role in host antimicrobial defense systems and stress responses 61 . Recently, it was regarded as a multi-functional protein, positioning MIF as a mediator during the inflammatory response to combat infection and in immunoinflammatory and autoimmune diseases 62 . At present, MIF homologs have been reported and studied in many crustacean species, such as S. paramamosain 63 , E. sinensis 64 , and Penaeus monodon 65 . In this study, MIF was significantly up-regulated at the post-molt stage, which may indicate that phagocytosis was promoted in S. paramamosain.
Cathepsin L (CTSL), a lysosomal cysteine protease, serves as a chemical barrier against microbial invasion in immune responses in vertebrates 66,67 . Cathepsin L has been found in crustaceans, such as E. sinensis 69 and F. chinensis 69 , is responsible for lysing pathogenic bacteria and as a barrier against invading pathogens.
C-type lectins (CTLs) are a large superfamily of pattern-recognition receptors that play important roles in the immune system through identifying and binding to the conservative pathogen-associated molecular patterns (PAMPs) on pathogen surfaces 70 . As pattern recognition receptors, CTLs can identify and bind PAMPs on pathogen surfaces to promote the phagocytosis 71,72 . In bivalve and crustacean species, the hepatopancreas synthesizes lectins (such as CTLs and ficolins) in response to infection 11,73,74 .
Hemocyanin, was reported as a novel and important defense molecule of the non-specific innate immune system [75][76][77] . Studies have found that hemocyanin functions not only as an oxygen-carrying protein, but also as a phenoloxidase-like enzyme 78 , and has antiviral activity, through its enzymatic cleavage, produces AMPs against a variety of viruses [79][80][81] . Agglutination activities against bacteria have also been observed as an antimicrobial strategy 82,83 .
To our knowledge, this is the first transcriptome analysis of the change of immune related genes in the hepatopancreas during the molt cycle in S. paramamosain (Fig. 6). Although these genes are usually multifunctional, they are significantly up-regulated at the post-molt stage and inter-molt stage revealing a potential role of the hepatopancreas in protecting crabs susceptible to infection. www.nature.com/scientificreports/  Sample collection. Mud crabs (S. paramamosain), averaging 7.9 ± 0.4 cm in carapace width and 100.2 ± 8.1 g in body weight, were collected from the aquafarm in Xiamen city, Fujian Province, China, and then transported to laboratory. They were reared in tanks (temperature: 26 ± 0.5 °C), fed with shrimp meat, and the water was changed every day. Mud crabs at three molt stages were sampled according to the previous study 16 . Post-molt stage: the mud crab just molted for 1-2 days, it was very thin, only dorsal carapace edge hard,Inter-molt stage: the dominant phase of the inter-molt period, the mud crab was full of meat, the exoskeleton was consolidated and hardened,Pre-molt stage: the fourth phase of the pre-molt period, namely double shell crab, pleural suture cracks between the carapace and the abdomen. Three individual crabs at each molt stage were sampled, placed on ice for anesthetization and then dissected for the hepatopancreas for RNA isolation. Three biological replicates were performed per tissue at each molt stage, for a total of 9 samples.
RNA isolation and cDNA library preparation. Total RNA was extracted from each tissue using TRIzol RNA isolation reagent (Invitrogen) according to the manufacturer's instructions. The RNA quantity and integrity were determined with an Aglient 2100 Bioanalyzer and a Qubit 2.0 Fluorometer. Magnetic beads with oligo (dT) were used to separate Poly-(A)-containing mRNA. Fragmentation buffer was added to randomly interrupt the mRNA into short fragments. The mRNA fragments were used as templates to synthesize the first-strand cDNA with random hexamer-primers, then add buffer, dNTPs, RNase H, and DNA polymerase I to synthesize the second-stranded cDNA, and further purified using Agencourt AMPure XP beads (Beckman Coulter). Next, double-stranded cDNA fragment ends were repaired, adding of poly (A) and ligation of adapters. The fragment size was selected using Agencourt AMPure XP beads. The PCR amplification to create the final cDNA library, and sequenced using the Illumina HiSeq2500.
cleaning and de novo assembly of sequencing reads. The assembles to compiled data from all 9 samples were assembled together for assembly, as this method is often more accurate and comprehensive for samples from the same species without a genome-wide reference 84 . The Trinity software first breaks the sequencing Reads into shorter fragments (K-mer), then extends the small fragments into longer fragments (Contig), and uses the overlap between the fragments to get the fragment. Finally, using the De Bruijn map method and sequencing Read information, the transcript sequences are separately identified in each fragment set 85

Differential expression analysis. Pearson's Correlation Coefficient (r) serves as an evaluation index
of biological repeat correlation 89 . DESeq 90 was used to identify the differentially expressed genes (DEGs). To ensure high quality of DEGs, Benjamini-Hochberg-corrected P values (q values) were calculated for one multiple hypothesis testing false discovery rate (FDR < 0.01) and fold change (FC ≤ 2) adopted as the key indicator for DEGs. The expression levels of up-regulated or down-regulated DEGs were annotated in the GO annotation database were enrichment analyzed using the topGO software 91 . COG database can perform orthologous classification of gene products and KEGG database helps to further interpret the function of gene. The Enrichment Factor was used to analyze the enrichment of pathway, and the Fisher's exact test method was used to calculate the enrichment significance.

Gene expression validation.
To validate and quantify the transcriptome data, genes identified in the DEGs were selected and then quantified using quantitative real-time PCR (qRT-PCR). Primers were designed using Primer 5.0 Tool (Premier, Canada). The stable housekeeping gene β-actin was used as an internal control (Table 4). qRT-PCR using a 7,500 Fast Real-Time PCR system (Applied Biosystems) was performed with the following thermal profile: One cycle at 95 °C for 30 s, followed by 40 cycles of 95 °C for 5 s, 58 °C for 30 s, and 72 °C for 30 s. The qRT-PCR was carried out in triplicate for each sample, and five mud crabs were analyzed in each group. Expression levels of different genes were estimated using the 2 -ΔΔCt method. Statistical significance (P < 0.05) was determined using Duncan's multiple range tests and one-way ANOVA under SPSS 18.0.