Transcriptomic analyses on muscle tissues of Litopenaeus vannamei provide the first profile insight into the response to low temperature stress

The Pacific white shrimp (Litopenaeus vannamei) is an important cultured crustacean species worldwide. However, little is known about the molecular mechanism of this species involved in the response to cold stress. In this study, four separate RNA-Seq libraries of L. vannamei were generated from 13°C stress and control temperature. Total 29,662 of Unigenes and overall of 19,619 annotated genes were obtained. Three comparisons were carried out among the four libraries, in which 72 of the top 20% of differentially-expressed genes were obtained, 15 GO and 5 KEGG temperature-sensitive pathways were fished out. Catalytic activity (GO: 0003824) and Metabolic pathways (ko01100) were the most annotated GO and KEGG pathways in response to cold stress, respectively. In addition, Calcium, MAPK cascade, Transcription factor and Serine/threonine-protein kinase signal pathway were picked out and clustered. Serine/threonine-protein kinase signal pathway might play more important roles in cold adaptation, while other three signal pathway were not widely transcribed. Our results had summarized the differentially-expressed genes and suggested the major important signaling pathways and related genes. These findings provide the first profile insight into the molecular basis of L. vannamei response to cold stress.


Introduction
The Pacific white shrimp, Litopenaeus vannamei, is one of the most important cultural shrimp species, and its production (by weight) has reached nearly 71% of the total economic penaeid shrimp production worldwide [1,2,3,4]. Naturally, L. vannamei can inhabit estuaries, lagoons or marine areas, but they are sensitive to low temperatures [1]. With high economic value, the farming areas of L. vannamei have been exploded from the southern to the northern coast, and a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 to the inland water in China. Therefore, low temperatures have become one of the major constrained factors to t L. vannamei culture, since that the culturing water temperature of the Northern China in winter is usually lower than that needed for L. vannamei survival.
Water temperature is one of the most common abiotic stress factors for aquatic ectotherms [5]. It influences nearly all biochemical and physiological progresses of aquatic organisms [6], and low temperatures may cause severe growth inhibition and mortality [6,7]. In fishes, the understanding of biochemical and molecular mechanisms under temperature stress have been well developed. In species such as rainbow trout (Oncorhynchus mykiss) [8], channel catfish (Ictalurus punctatus) [9] and common carp (Cyprinus carpio L.) [10], studies showed that about 10 percent of the investigated genes were found to be up-regulated under temperature stress, and they were involved in an extensive gene expression network [11]. Digital gene expression in Nile tilapia (Oreochromis niloticus) revealed the differentially expressed-genes when comparing the fish cultured in 10˚C and 30˚C [12]. These studies provided useful insights and promoted the molecular breeding work in fish species. However, in L. vannamei, the gene expression profile in response to temperature stress has not yet been illustrated.
RNA-Seq technology is a crucial method for quantifying the transcriptional expression in non-model species or in those organisms that lack whole-genome sequencing data [13,14]. Studies on channel catfish (I. punctatus) have suggested that both the de novo synthesis of cold-induced proteins and the modification of existing proteins were required for adaptation to low temperatures [9]. Transcriptional regulation research in the coral reef fish Pomacentrus moluccensis revealed the existence of differentially expressed genes in this species [15]. Moreover, RNA-Seq has already been applied as a useful method on other stress response studies in L. vannamei, such as high salinity [16] and low salinity [17,18,19,20]. However, the specific expression profile report on L. vannamei under the stress of cold temperature is still lacked.
Molecular signaling pathways can help to improve the understanding of the activity and coordination of the cell-expression network in response to extracellular stress. Calcium is a ubiquitous second messenger in eukaryotic signal transduction cascades [21,22]. It impacts nearly every aspect of cellular life [23] and fulfills its key role in transferring extracellular stress signals into chromatin and in regulating protein phosphorylation [24]. Structural comparisons of mammalian Ca 2+ /calmodulin-dependent protein kinase II (CaMKII) and plant calciumdependent protein kinases (CDPKs) of calcium signaling showed that their kinase domain containing the highly conserved subdomains of typical eukaryotic Ser/Thr protein kinases [22]. Mitogen-activated protein kinases (MAPKs), regulated by MAPK kinase kinase-MAPK kinase-MAPK (MKKK-MKK-MAPK) phosphorelay system, are Ser/Thr protein kinases and are key components for vital signal transduction pathways that converting extracellular stimuli into a wide range of cellular responses in the organisms from yeast to humans [25,26,27]. Consequently, in order to understand the molecular basis of L. vannamei under cold temperatures, investigations conducted on differentially expressed genes in signaling were necessary.
Shrimps are poikilothermic species, for which changing of water temperatures can conduct to each part of their bodies and influence their cellular activities. Muscle compose nearly 60% body weight of L. vannamei, and it is the main tissue for converting biochemical energy in shrimp [6]. Furthermore, we found that the muscle tissues of L. vannamei turned whitish immediately when they were chilled in low temperature water, but regained pellucidity after being returned to the control temperature, indicating that the muscle tissues of L. vannamei may process overt physiology activities in response to acute cold stress. In the present study, a transcriptomic analysis has been carried out on the muscle of L. vannamei. Gene expression profiles under the stress of low temperature has been characterized by comparing various libraries from corresponding temperature treatments, including acute cold stress for 2 hours, endurance cold stress for 48 hours and recovering to normal temperature for 2 hours. The identified candidate genes will be useful and helpful in revealing a deeper insight into the molecular adaptation basis under low temperature stress in L. vannamei.

Ethics statement
Shrimp care and experiments were carried out according to the Care and Use of Agricultural Animals in Agricultural Research and Teaching and approved by the Science and Technology Bureau of China. Approval from the Department of Wildlife Administration was not required for the experiments conducted in this paper. All experiments in this paper were performed with permits obtained from the Government of the People's Republic of China and endorsed by the Animal Experimentation Ethics Committee of Chinese Academy of Sciences.

Treatments and cold challenges of the shrimp materials
Healthy shrimp with an average weight of 7.09±3.22 g were treated in gradient temperatures (9,10,11,12,13,14,15,18,21,24,27,30 and 33˚C) for 2 hours, 48 hours and then recovered at a control temperature (28˚C) for 2 hours. Survival was evaluated by observing the shrimp's ability to move spontaneously or after gentle prodding [28]. The temperature of semi-death group was set as the stress cold temperature.
Thirteen degree was set as the stress cold temperature. Four groups were placed either in control temperature (CT), 13˚C for 2 hours (T_2), 13˚C for 48 hours (T_48) and then, the T_48 group was placed to recover temperature for 2 hours (RC). Muscle tissues from the surviving individuals were put in centrifuge tubes and dipped into liquid nitrogen immediately. Three biological replicates in equal amounts of each group (2 μg of RNA for each individual, and total 6 μg of RNA for each group) were pooled for subsequent RNA-Seq experiments.

RNA isolation and cDNA library construction
Four separate libraries (CT, T_2, T_48 and RC) were sampled. Total RNA was extracted by TRIzol Reagent (Invitrogen, USA) according to the manufacturer's instructions. The quality of the RNA products was verified as described by Ren et al. (2014) [14]; RNase-free DNase I (Invitrogen, USA) was used to degrade any possible DNA, and oligo (dT)-coated magnetic beads were mixed with the total RNA to concentrate the ployA tailed mRNA. Fragmentation was performed by incubation in an NEB Next First Strand Synthesis Reaction Buffer (NEB, USA), and the second strand was generated with the buffer, dNTPs, RNase H and DNA polymerase-I. Adapters were ligated to the synthesized cDNA fragments after an end repair step.

RNA-Seq and annotation
The library was constructed using an Illumina HiSeq 4000 sequencing platform located at the BGI Company (Shenzhen, China; http://www.genomics.cn/index). The clean read data were deposited to the US National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA, http://www.ncbi.nlm.nih.gov/Traces/sra) with the accession No. SRP095377.

Quantification of differentially expressed genes
Three comparisons were taken to the four libraries, including CT-VS-T_2 (T_2/CT), T_2-VS-T_48 (T_48/T_2), and T_48-VS-RC (RC/T_48). The differentially expressed genes (DEGs) in comparison of T_2/CT represented the transcripts in response to acute cold stress (placed from control temperature to acute 13˚C). The DEGs in comparison of T_48/T_2 represented the transcripts in response to endurance cold stress (endured the cold stress for 48 hours). The DEGs in comparison of RC/T_48 represented the transcripts for recover adaptation (placed from stress cold to control temperature). The clean reads per kilobase per million (RPKM) value was used to estimate transcript abundance on the basis of eliminating the influence of different gene lengths and sequencing discrepancies [31]. Therefore, the RPKM values could be directly used to compare the differences of gene expression among samples. The algorithm developed by Audic and Claverie (1997) [32] was used to compare the differences in various gene expression levels. A value of the FDR (false discovery rate) of less than 0.001 and a value of the Log2 ratio of less than 1 were set as the criteria values.

Identification of major response genes and pathways
The coexisted DEGs in different comparisons of total genes, GO annotated genes and KEGG annotated genes were filtered. Four important kinds of signaling genes regarding Calcium and MAPK cascades, Transcription factor and Serine/threonine-protein kinase signaling genes were searched from the total annotated genes [14]. To understand the expression model of target genes, hierarchical clustering analysis was performed with the software Heatmap Illustrator Heml 1.0.3.3.

Real-time PCR validation
To verify the expression profiles of genes in our RNA-Seq results, the represented 15 DEGs in comparison of T_2/CT were selected for validation of the Illumina sequences by real-time PCR analysis. All primers used were listed in Table 1. The treatments of L. vannamei, sampling of muscle, isolation of total RNA and treatment of DNase were performed as described above. The first cDNA strand was synthesized by PrimeScript1 Reverse Transcriptase Kit (Takara, China) with an oligo (dT) primer, a qPCR analysis was conducted in a RotorGene RG3000 real-time PCR system (Corbett research, Australia), PCR reactions were conducted using a SYBR Premix Ex TaqTM II (Takara, China), and the reactions were exposed to an initial denaturation (94˚C for 3 min) followed by 40 cycles of 94˚C for 25 s, 60˚C for 15 s, and 72˚C for 15 s. The relative transcript abundance was obtained by normalizing with expression of the L. vannamei β-actin gene based on the 2 -ΔΔCT method [33].

Cold stress temperature
A batch of L. vannamei were treated with different temperatures, and the survival rates were summarized ( Fig 1B). The results showed that L. vannamei could not survive when treated with 11˚C or lower temperatures, the survival rates were slightly affected in the water at 14˚C, 15˚C or higher temperatures ( Fig 1B). The survival rates at 12˚C and 13˚C were 9.30% and 47.22%, respectively. We used 13˚C as the cold stress temperature since a nearly half survival rate (47.22%) of shrimp was observed at this temperature. The same batch of shrimp from the same cultural environment was used for the following experiments.

Overview of the libraries
After filtering, the quality metrics of the reads are showed in Table 2: the total raw reads of the samples were from 87.69 to 90.94 Mb, the number of clean reads per library ranged from 65.23 to 66.02 Mb, and the clean reads ratio ranged from 71.95% to 75.28%. The mean value of clean reads Q20 and Q30 percentages were 97.41% and 94.34%, respectively. The clean reads were assembled and clustered into Unigenes, and the quality metrics of Unigenes are presented in Table 3. The total number of All-Unigenes was 29,662, the total length was 29,632,338 bp and the mean length was 999 bp. The annotated genes in NR, NT, SwissProt, KEGG, COG, GO database were 17,103 (57.66%), 11,976 (40.37%), 15,552 (52.43%), 13,936 (46.98%), 8,558 (28.85%) and 2,504 (8.44%), respectively (Fig 2A). Overall 19,619 Unigenes were annotated (Fig 2A). The main Venn diagram of NR, KEGG, COG and SwissProt annotation are showed in Fig 2B, the largest interior gene numbers were 8,264. The percentage of annotated species with NR database are displayed in Fig 2C, the top three annotated species were Zootermopsis nevadensis (the largest blue part, 12.34%), Daphnia pulex (the second large tenne part, 5.31%) and Tribolium castaneum (the third large gray part, 2.95%), respectively. GO and KEGG classification of the DEGs GO annotation was conducted for the three comparisons. In the comparison group of T_2/ CT, 870 of 1471 DEGs could be assigned by GO classification, and the equivalent numbers for other comparison groups were 1147 of 1943 DEGs (T_48/T_2) and 690 of 1477 DEGs (RC/ T_48); the detailed GO classification data (Pathway, PathwayID, and GeneID) are listed in S1 Table. The most annotated GO pathway in each comparison was Catalytic activity (GO: 0003824). KEGG annotation was conducted for the three comparisons. In the comparison of T_2/CT, 1202 of 1471 DEGs were assigned to 83 pathways. In the comparison of T_48/T_2, 1556 of  1943 DEGs were assigned to 91 pathways. In the comparison of RC/T_48, 1208 of 1477 DEGs were assigned to 77 pathways. The details of the KEGG classification (Pathway, PathwayID, and GeneID) of the three comparisons are listed in S2 Table. The top two annotated KEGG pathways in all the three comparisons were Metabolic pathways (ko01100) and RNA transport pathways (ko03013). The RNA transport pathway maps (there were no Metabolic pathway maps in KEGG database) were traced and downloaded from the KEGG database (S1 Fig), the annotated DEGs in the three comparisons were summarized (S3 Table), the major discrepant DEGs in RNA transport pathway were the components of exon-junction complex (EJC). The annotated up-and down-regulated genes of EJC complex were marked and displayed in Fig 3. The EJC complex mainly participated in mRNA translation (S1 Fig, Fig 3A), the EJC outer shell and Transiently interacting factor genes were the most two kinds of responding components in EJC complex (S3 Table,

Quantification of transcripts
The transcript abundance of each Unigene from the four related libraries is listed in S4 Table. DEGs were identified through a pairwise comparison between various transcripts by setting a threshold FDR (false discovery rate) value of 0.001 and a |log 2 The down-regulated transcripts were far more than the up-regulated genes in comparison of T_48/T_2, indicating that more genes were suppressed under the endurance cold stress.

Identification of major response genes
The Venn diagrams in Fig 6A were used to exhibit the coexisted DEGs in different comparisons. It showed that 234 genes were coexisting in the three comparisons of Total DEGs, and 35 of coexisted GO annotated DEGs and 48 of coexisted KEGG annotated DEGs were also identified.
Top 20% of differently-expressed coexisted genes in the three comparisons were selected and taken to Hierarchical clustering analysis, 72 genes were obtained, the detailed GeneID and other information are showed in S8 Table,

Identification of major pathways
Total 35 coexisted GO-annotated DEGs and 48 coexisted KEGG-annotated DEGs were selected, the detailed GeneID information are showed in S9 and S10 Tables. The major (Absolute Value of log 2 FoldChange ! 5.0) Unigene in GO and KEGG Pathways were filtered, and 15 GO pathways and 5 KEGG pathways were recognized as the temperature-sensitive pathways in response to cold stress ( Table 4). The Unigene of CL596.Contig2_All, described as Eukaryotic translation initiation factor 5 gene (gi|646705355|), was annotated by Translation initiation factor 5 pathway (KEGG: K03262) and Cellular macromolecule metabolic (GO: 0044260) and primary metabolic process (GO: 0044238), implying that this gene might play important role in several molecular functions in response to cold stress.

Major signaling genes in response to cold stress
Four important kinds of signaling genes regarding Calcium signaling, MAPK cascades, Transcription factors and Serine/threonine-protein kinase were fished out (S11-S14 Tables) and clustered into heat map. Generally, the Calcium signaling genes, MAPK cascades genes and Transcription factor genes were not widely transcribed in response to cold stress (Fig 6C, 6D and 6E). Four Unigenes of Calcium signaling genes, 4 Unigenes of MAPK cascades genes and 5 Unigenes of Transcription factor genes were obtained when filtered by the top two levels of log 2 FoldChange (absolute value of log 2 FoldChange ! 6.0 in Calcium signaling genes, absolute value of log 2 FoldChange ! 5.0 in MAPK cascades genes and absolute value of log 2 FoldChange ! 3.5 in Transcription factor genes, Table 5). Whereas, conspicuous patches of Serine/threonine-protein kinase genes were observed (Fig 6F), 16 Unigenes were obtained when filtered by the top two levels of log 2 FoldChange (absolute value of log 2 FoldChange ! 5.0). The results suggested that Serine/threonine-protein kinase might participate in the adaptation of acute or endurance cold stress.

Real-time PCR validation
The  qPCR and RNA-Seq were also displayed, and the results showed all the 15 candidate genes in qPCR verification consisted with the results of the RNA-Seq technology (Fig 7B).

Discussion
Temperature is one of the most important environmental variances for aquatic ectotherms, and it has been investigated for its key role in affecting overall biological processes in various species [7,34]. In larval zebrafish (Danio rerio), cold stress induced-transcripts were generated from 2 and 48 hours treatments under 16˚C [11]. In tilapia (O. niloticus), an analysis of gene expressions was carried out under cold stress by comparing different temperatures, from 30˚C followed by a decrement of 1˚C/Day to 10˚C [12]. In Manila clam (Ruditapes philippinarum), a transcriptome study was conducted in response to cold stress with an extreme low temperature of -1˚C [35]. Stresses could be classified into mild, chronic or acute, which represented a dramatic shift in the environmental conditions [36]. In our study, shrimp were treated in acute cold temperatures for 2 hours, 48 hours and then in a recovery temperature for 2 hours, and a nearly 50% survival rate (47.22%) was observed in the 13˚C group. The stress temperature in our result was similar to that reported by Fan et al. (2013) [37], in which 13˚C for 24 hours was regarded as a suitable cold stress temperature to investigate the altered proteins and genes in L. vannamei.
RNA-Seq technology has been demonstrated as a straightforward method to reveal the general molecular regulation processes in species underlying the response to stressful environments [14]. In the previous studies for L. vannamei, RNA-Seq technology has been employed to investigate the different gene expression levels involved in various environment factors and developmental processes, such as a transcriptome analysis under acute ammonia stress [38], some digital gene expression analyses in response to low salinity stress [17,18,19,20,39], and the transcriptional profiling during ontogenesis [40]. Parts of the molecular mechanisms have already been characterized in response to temperature in the shrimp. The protein content of CAT, GST, Ferritin and HSP60 have been found to be influenced by temperature stress in L. vannamei [41]. Lipid peroxidation, catalase activity and glutathione-S-transferase in Palaemon elegans and Palaemon serratus have been proved as the biomarkers for monitoring aquatic environments [42]. In addition, the decreasing water temperatures has been shown to induce oxidative stress in L. vannamei [43,44]. Furthermore, 20 proteomic spots and several genes have demonstrated to play important roles in temperature regulation of L. vannamei [37]. However, the profile information of the regulation of gene expression under low temperatures in L. vannamei was still lacked. In the present study, transcriptomic analyses were conducted under cold temperature stress on the muscle of L. vannamei. To our knowledge, this is the first report of transcriptional profiling under acute low temperature stress in L. vannamei.
Transcriptome analyses may provide insights into the molecular basis of cold tolerance in many species. Transcriptomic characterization under temperature stress in zebrafish (D. rerio) revealed the transcriptional and post-transcriptional regulation of the cold-acclimated processes [5], the multiple temperature-regulated biological processes and pathways [11], and the over-expressed biological processes (such as circadian rhythm, energy metabolism, lipid transport and metabolism) [19]. An RNA-Seq analysis in Melanotaenia duboulayi provided candidate genes in response to changing temperatures in freshwater fish [45]. The study of transcriptomic analysis on the marine fish ectoparasitic ciliate Cryptocaryon irritans under a low temperature revealed the required genes for cell survival and a deeper dormancy state [46]. In our present study, a transcriptomic analysis for L. vannamei under acute low temperature was carried out. In total of the four libraries, 904 up-and 567 down-regulated genes in comparison of T_2/CT, 634 up-and 1309 down-regulated genes in comparisons of T_48/T_2, and 997 up-and 480 down-regulated genes in comparisons of RC/T_48 were detected. The top 20% of coexisted DEGs were fished out. In our result, 15 GO pathways and 5 KEGG pathways were annotated as temperature-sensitive pathways. Our study promoted the understanding of profiling gene expression of L. vannamei in response to cold stress.
Acute stress places cells in danger, and rapid adaptation is crucial for maximizing cell survival [36]. To adapt to either biotic or abiotic stress, cellular signaling pathways may coordinate multiple integrated and regulated layers of different steps of mRNA biogenesis in inner cells [24,36]. Thus, the important signaling pathways regarding cellular stress response adaptation have been investigated in many studies. Calcium ions (Ca 2+ ) impact nearly every aspect of cellular life [23]; Teets et al. (2013) demonstrated that goldenrod gall fly (Eurosta solidaginis) tissues used calcium signaling to instantly detect decreases in temperature and to trigger downstream cold-hardening mechanisms [47]. In the transcriptomic analysis reported by Ren et al. (2014), the gene members of the calcium signaling, including cation/calcium exchangers, calcium-binding proteins, calmodulin-like proteins, etc., were found to be involved in the response to cold stress in Chrysanthemum morifolium [14]. The MAPKs influenced the cellular signal transduction in organisms from yeast to human. When yeast under osmostress, the Ste11 MAP3Ks were activated, followed by activation of Pbs2 MAP2K that combined to Hog1 and then phosphorylation of the specific osmostress transcription factors for cell adaption [36]. The p38 family in mammalian species contains 4 isoforms (MAPK11-14) with many overlapping functions in cells, and downstream targets of p38 MAPKs include several kinases that are involved in the control of gene expression and nuclear proteins, such as transcription factors (CREB, ATF1, NFκB, STAT1 and STAT3) and regulators of chromatin remodeling (H3 and HMG14) [36]. In this paper, the expression model for the genes in four signal pathways were investigated under cold stress. The major differentially-expressed signaling genes were suggested involve in the adaptation of cold stress. Serine/threonine-protein kinase genes might play more important roles in response to cold stress, while Calcium, MAPK cascades and Transcription factor signaling genes were not widely transcribed in response to cold stress.
In conclusion, the overall transcriptional expression and quantification analyses under the cold stress temperature of 13˚C were conducted for Litopenaeus vannamei. Top 20% of differentially expressed genes were filtered. The most annotated and temperature-sensitive of GO and KEGG pathways were found out. The major differentially-expressed genes in signaling were suggested. The results showed that Ser/Thr kinase signal pathway might play more important roles in participating the cold adaptation, while Calcium, MAPK cascades and Transcription factors signaling genes were not widely transcribed. Our study provided the first insight into the molecular basis and supplied an important reference for the mechanism of L. vannamei responded to cold stress.  Table. Different transcribed genes between the comparison of T_48/T_2. a The "Gene-Length" column gives the length of exon sequence. b T_2: reads per kb per million reads (RPKM) for each unigene in transcription of T_2. c T_48: reads per kb per million reads (RPKM) for each unigene in transcription of T_48. d Log2FoldChange(T_48/T_2): the ratio between the RPKM in T_2 to the RPKM in T_48. The criteria applied for assigning significance were: Pvalue < 0.05, FDR 0.001, and estimated absolute |log2FoldChange| ! 1. Genes listed in descending order of absolute |log2FoldChange|. NA: No annotation were found in the related database.