Analysis of Stress-Responsive Transcriptome in the Intestine of Asian Seabass (Lates calcarifer) using RNA-Seq

Identification of differentially expressed genes (DEGs) and regulated pathways in response to stressors using a whole-genome approach is critical to understanding the mechanisms underlying stress responses. We challenged Asian seabass with lipopolysaccharide (LPS), Vibrio harveyi, high salinity and fasting, and sequenced six cDNA libraries of intestine samples using Roche 454 RNA-seq. Over 1 million reads (average size: 516 bp) were obtained. The de novo assembly obtained 83 911 unisequences with an average length of 747 bp. In total, 62.3% of the unisequences were annotated. We observed overall similar expression profiles among different challenges, while a number of DEGs and regulated pathways were identified under specific challenges. More than 1000 DEGs and over 200 regulated pathways for each stressor were identified. Thirty-seven genes were differentially expressed in response to all challenges. Our data suggest that there is a global coordination and fine-tuning of gene regulation during different challenges. In addition, we detected dramatic immune responses in intestines under different stressors. This study is the first step towards the comprehensive understanding of the mechanisms underlying stress responses and supplies significant transcriptome resources for studying biological questions in non-model fish species.


Introduction
Stress response is important to living organisms and a major aspect of natural selection in the wild. 1 Fishes are exposed to a variety of stressors because their homeostatic mechanisms are highly dependent on prevailing conditions primarily related to the aquatic environment. 2 Previous studies showed that the interactions between stressors and stress responses are complex. 2 Genetic, developmental and environmental factors, and early experiences have influences on the stress responses of fish. 3,4 Identification of differentially expressed genes (DEGs) and pathways under stress conditions is important for the understanding of fish stress response. A lot of platforms to reveal DEGs are available, such as real-time (RT)-PCR and microarray. Recently, RNA-Seq 5 is considered to be a revolutionary tool for transcriptomics, as it can absolutely quantify over millions of unknown transcripts. It has also shown a great analytical power in the identification of differentially expressed transcripts in response to different conditions. 6 -9 For species with poor genome annotation, the Roche/454-based RNA-seq technique with much longer read lengths provides an excellent tool for expression profiling studies. 10,11 The intestine is one of the major organs in fish that interact with the environment and is involved in many biological processes. A few studies on the responses of the intestine to various stressors have been performed in fishes, such as medaka 12 and channel catfish. 13 However, data on responses of the fish intestine to different stresses are still limited. Studies of the genome-wide transcriptome after challenging with different stressors in fish will enhance our understanding of the molecular mechanisms underlying stresses.
The Asian seabass (Lates calcarifer) is an important food fish species in Southeast Asia and Australia. 14 Although some expressed sequence tags (ESTs) of the Asian seabass are available in public domains, 15 -17 the transcriptome for the Asian seabass is still poorly characterized. To date, an intestine transcriptome data set for the Asian seabass and most food fish species is lacking. To gain more understanding of the mechanisms underlying stress responses and to supply significant transcriptome resources for studying biological questions in non-model fish species, we conducted RNA-seq of intestine samples from fishes challenged with lipopolysaccharide (LPS), Vibrio harveyi infection, fasting and high salinity. We also carried out bioinformatic analyses of the transcriptome to identify DEGs and pathways in response to these different stressors.

Fish management, challenges and sampling
for RNA-seq analysis Thirty-six Asian seabass at the age of 11 months (body weight 330 g) were originally maintained in a large tank containing 2000 l of freshwater in the animal outhouse of our institute. For challenge experiments, 12 fishes were transferred to a 1000-l tank, and the salinity concentration was gradually increased to full-seawater (33 PPT salinity) within 3 days. Fishes were fed twice daily with pelleted feed (Biomar, Nersac, France). One day prior to challenges, nine seabass from the seawater tank, after acclimatization for 2 weeks, were evenly divided into three tanks containing 300 l of seawater (i.e. 3 fishes per tank). For the three fishes in tank 1, named as Int1 (LPS), each fish was injected intra-peritoneally with 0.3 ml of 5 mg/ml of Escherichia coli LPS (Sigma-Aldrich, Saint Louis, USA) by dilution with phosphate-buffered saline (PBS) at room temperature. In tank 2, Int2 (Vibrio), a total of three fishes were injected intraperitoneally with 0.3 ml of PBS dissolved culture pellet of Vibrio harveyi ( e 8 cell/ml) at room temperature. In tank 3, Int3 (PBS), used as control 1, three fishes received an intra-peritoneal injection of 0.3 ml of PBS for each fish. These fishes were not given access to feeds before sampling. Three fishes taken from the original freshwater tank were moved to tank 4, Int4 (FW;Feed), as control 2, containing 1000 l of freshwater. These fishes were fed twice daily with pelleted feed (Biomar, Nersac, France). Another three fishes from the original freshwater tank were moved to the freshwater tank 5, Int5 (FW;Fasting), and were not given access to feed before sampling. Three fishes in the seawater tank 6, Int6 (SW;Feed), were fed twice daily with pelleted feed before sampling. Three fishes from each of the tanks 1, 2 and 3 were sacrificed at 40 h post-challenges. Three fishes from each of the tanks 4, 5 and 6 were sacrificed at 8 days posttreatments. Intestine samples were taken from each fish of each tank and kept in Trizol reagent (Invitrogen, Carlsbad, USA) for RNA isolation.

Challenges and sampling for quantitative RT-PCR
analysis Eighteen seabass at the age of 11 months (body weight 330 g) originally maintained in a large tank containing 2000 l of freshwater were evenly divided into two tanks containing 1000 l of freshwater (Groups 1 and 2). Nine of the fishes in Group 1 were not given access to feed before sampling, and nine of the remaining fishes in Group 2 were fed twice daily with pelleted feed (Biomar, Nersac, France). Three fishes from each group were sacrificed at 8 days postfasting. Intestine samples were taken for each fish and kept in Trizol reagent (Invitrogen, Carlsbad, USA) for RNA isolation.
For analysis of the functions of the splice variants of IFABP-a and -b genes, two extreme groups (i.e. smallest and biggest; n ¼ 6/group) for body weight were selected from a population of 300 seabass at the age of 2 months. These fishes were originally maintained in a tank containing 2000 l of freshwater and were fed twice daily with pelleted feed (Biomar, Nersac, France). Intestine samples were taken for each fish and kept in Trizol reagent (Invitrogen, Carlsbad, USA) for RNA isolation.

RNA-seq sequencing
Total RNA from the intestine was isolated using the Trizol kit (Invitrogen, Carlsbad, USA). Total RNA from three fishes at each time point was equally mixed and submitted to Macrogen (Seoul, Korea) for RNA sequencing by using Roche/454 GS FLX Titanium platform. The total RNA quality was assessed using the Agilent 2100 Bioanalyzer. Ribosomal RNA was then removed prior to proceeding. cDNA rapid libraries were prepared according to the manufacturer's protocol (Roche, Central plaza, Singapore).

2.4.
De novo assembly of the intestine transcriptome for the Asian seabass GS FLX data were processed using the Roche GS FLX software (v 2.6). De novo assembly of transcriptome was carried out using the GS De Novo Assembler Singletons and isotigs from the de novo assembly were referred to as unique sequences to annotate them. They were blasted against GO database (http://www. geneontology.org/), Swissprot, Ref_protein and Refseq_RNA databases that retrieved from NCBI database using Blastn or Blastx algorithm with an E-value threshold of E 25 .

Annotation of pathways
Ortholog assignment and mapping of the unisequences to biological pathways were performed using the KEGG automatic annotation server (KAAS; http://www.genome.jp/ tools/kaas/) 18 with a threshold bit-score value of 40.

Prediction of protein-coding regions and se-
quence alignment OrfPredictor 19 was used to predict protein-coding regions in EST-derived unisequences. The predicted protein sequences were then input into SignalP 4.0 server 20 to predict potential signal peptide cleavage sites. Alignment of the predicted protein sequence for each splice variant was conducted using the online tool ClustalW2 (http://www. ebi.ac.uk/Tools/msa/clustalw2/).

Identification of DEGs and significantly regu-
lated pathways Identification of DEGs from the count data was performed using the NOISeq-sim program 21 with parameters setting as v ¼ 0.1, k ¼ 1.5, and other as defaults. Normalization was conducted using the RPKM method. 22 The identified DEGs were then blasted against the KAAS database to identify the regulated pathways involved. Comparison of the DEGs and pathways among the four challenges was carried out with Venn diagrams. 23 The RPKM-normalized transcript counts for the top 60 most DEGs in each pairwise comparison were then converted to natural logs and analysed by Cluster 3.0 with parameters as hierarchical clustering, uncentered correlation, and complete linkage (http://bonsai.ims.u-tokyo.ac.jp/~mdehoon/ software/cluster/software.htm#ctv) to construct a heat map. Scatter plots between pairwise data sets were drawn with geWorkbench platform (https:// cabig.nci.nih.gov/community/tools/geWorkbench).
2.5.5. qPCR analysis of gene expressions and functions of splice variants of IFABP genes For quantitative RT-PCR (q-PCR) analysis to confirm gene expression revealed by RNA-seq, total RNA from the same group (experimental group: three fishes sampled at 8 days post-fasting challenge and the control group) were equally mixed. Around 1 mg of DNase I-treated total RNA was reverse transcribed to cDNA by M-MLV reverse transcriptase (Promega, Madison, USA) with random hexamer primer as RT primer following the manufacturer's protocol. The resulting single-strand cDNA was assayed as DNA template by qPCR using primers for 56 genes from the RNA-seq data and the EF1A gene as the control (Supplementary Table S2). PCR in triplicates was performed with the KAPA TM SYBRw FAST qPCR Kits (Kapa Biosystems, Boston, USA) as described by the manufacturer, and an iQ5 PCR-machine (Bio-Rad, CA, USA). The DDCT method was used for the analysis of the gene expression, and the values of triplicate RT-PCR reactions were normalized to EF1A gene expression.
The qRT-PCR analysis of expressions of the splice variants of two IFABP genes in the intestines of the Asian seabass were performed as described above with slight modifications. Briefly, to design primers for PCR, the exon/intron boundaries and divergence regions from the splice variants of IFABP genes were identified by using Sequencher 4.9 (GeneCodes, CA, USA). One primer pair that spans exon/intron boundary on the mRNA for each splice variant was designed using PrimerSelect (DNASTAR, DE, USA). This primer design allows differentiating the splice variants and does not amplify their genomic DNA. In addition, around 1 mg of DNase I-treated total RNA from each fish was reverse transcribed to cDNA. The resulting singlestrand cDNA was assayed as DNA template by qPCR using primers for the six splice variants. The EF1A gene was used as a control (Supplementary Table S2).

Statistical analysis
The correlation coefficients between pairwise data sets were calculated by using the SAS (SAS Institute, Cary, USA). A P-value for significant difference was calculated using the t-test module installed in the Microsoft Office Excel 2008 program with parameters setting as two-tailed distribution and two-sample equal variance.

Results and discussion
3.1. De novo assembly and annotation of an intestine transcriptome for the Asian seabass 3.1.1. Transcriptome sequencing by using 454 technology and de novo assembly A total of 1 004 081 raw sequencing reads with an average length of 516 bp were generated by sequencing 6 cDNA libraries. The high-quality reads generated are available at the NCBI SRA browser (Run accession number: DRR002185-DRR002190). The de novo assembly produced a total of 83 911 unique sequences (referred to as unisequences) including 33 191 isotigs and 50 720 singletons ( Table 1). The average size of isotigs was 992 bp and the N50 isotig size was 1191 bp, while the average length of the unisequences was 747 bp (Supplementary Fig. S1a and b). Compared with other de novo transcriptomes in fish, such as spleen transcriptome of the orange-spotted grouper with an average length of 504 -547 bp, 24 our study produced a much better transcriptome assembly. Contigs and singletons are also available upon request.
By using Blastn search reciprocally, only 46.46% (38 992) of the unisequences in the new built had significant matches to the Asian seabass NCBI EST data set. The remaining 44 919 (53.54%) unisequences without significant similarity to the previous EST data set were considered to represent new sequences in Asian seabass. We also compared the new assembly with the unigene sequences of the model fish species to evaluate the depth and breadth of our unisequence data set. 40.82 and 67.41% of transcriptome data sets of Danio rerio and Gasterosteus aculeatus had significant hits in the seabass intestine transcriptome, respectively. However, only 6.63 -17.00% unisequences of the intestine transcriptome could be mapped onto the unigenes sequences of Takifugu rubripes and Gasterosteus aculeatus, respectively. These data suggest that a large number of new gene sequences are discovered in the Asian seabass intestine transcriptome. This intestine transcriptome can serve as a reference sequence database for intestine comparative transcriptome analysis in fish. The summary of reciprocal comparison between public EST data sets and the intestine transcriptome data set of the Asian seabass by Blastn program is presented in Supplementary Table S1.  Table S5). Because a signal peptide directs the transport of a protein, therefore, 4.48% of the transcriptome proteins were likely destined to the secretory pathway in the fish intestine. Our predictions are similar to the published data of other species, such as parasites (3305 -4246 sequences; 6.5 -6.9%) 25 and trematode (1534 sequences; 5.03%). 26 3.1.4 Verification of differential expressions of genes by qPCR To validate the digital gene expression analysis from the intestine RNA-seq data, 56 genes were randomly selected for qPCR analysis of samples at 8 days post-fasting challenge and the control. In the two data sets, most genes showed similar responses to the fasting challenge. The RNA-seq data had a positive linear relationship with qPCR data (Pearson correlation coefficients (r) ¼ 0.58; Supplementary Fig. S2). There was no statistically significant difference (t-test: P ¼ 0.21) between two data sets, which was in concordance with the report in channel catfish RNA-seq study. 13 These results indicate that the RNA-seq technique could provide an excellent tool for expression profiling studies. For non-model organisms without reference genome and transcriptome sequences, although RNA-seq with the Illumina sequencing platform could generate a large number of short sequence reads, it is very difficult to assemble and analyse the reads. Therefore, 454 sequencing generating much longer reads is the suitable choice in Asian seabass. It is to note that due to a lower number of reads in 454 sequencing in comparison with the Illumina sequencing, some extremely low expression genes might not have been detected in this study. In the future, RNA-seq with the Illumina sequencing technology could be used for more detailed analysis of profiles of gene expression.

Stress responses in the Asian seabass intestine
3.2.1. Identification of DEGs and regulated pathways in response to challenges The interactions between stressors and stress responses are highly complex, and some stress responses may themselves act as stressors and vice versa. 2 To obtain a general view of DEG expression patterns, pairwise comparison of expression abundance in the RNA-seq data sets was first conducted. The scatter plots showed a linear relationship between pairwise data sets with positive Pearson correlation coefficients (r ¼ 0.82 to 0.99, P , 0.01) ( Supplementary Fig. S3), suggesting that expression patterns of most genes are similar, and only a small portion of genes were expressed differentially among different stresses.
To explore DEGs in response to different challenges, we performed statistical analysis of gene expression between pairwise samples using the NOISeq procedure. Because the fishes used in this analysis are heterogeneous populations, more restrictive parameters than the defaults in the software were applied during data analysis. Volcano plots in Fig. 1 show gene expression differences between challenges and controls. The numbers of DEGs in response to LPS and Vibrio infection were higher than those induced by fasting and highsalinity stressors. These results suggest that some molecular response mechanisms are dependent on the nature of the stress signals, and LPS and Vibrio harveyi infection are the stronger stressors among the four different challenges (i.e. LPS, Vibrio harveyi, high salinity and fasting).
When using probability 0.8 (means that the feature is four times more likely to be differentially expressed than non-differentially expressed) as the threshold for each gene of being differentially expressed, we found 1080 (Fasting) to 1739 (LPS) genes were differentially expressed between treatments and the respective controls (Table 2). Of which, 601 (high salinity) to 1139 (LPS) DEGs were down-regulated, and 410 (Fasting) to 627 (high salinity) DEGs were up-regulated. The observed values for DEGs are similar to those in the previous study based on RNA-seq in channel catfish, in which, 1633 DEGs at 3, 24 h, and 3 days between challenged and control intestine samples were revealed following Edwardsiella ictaluri infection. 13 The DEGs in this study are similar to those in catfish. 13 For example, around half of the DEGs in the category Cytoskeletal/ Muscle Fiber Dynamics that were detected in channel catfish can be found in our DEG category following Vibrio infection, such as Actin (cytoplasmic 1), AHNAK nucleoprotein and Annexin A2. These observations suggest that some physiological responses to stressors in fish are well conserved evolutionarily. Around 31% (Vibrio infection) to 50% (LPS injection) of the DEG sequences in this study had no hit in the current protein or RNA databases, suggesting that these genes are novel. These newly identified DEGs are of interest for further study. The list of DEGs that were identified using NOISeq ( probability 0.8 as a threshold) is presented in Supplementary Table S6. To illustrate the differential expression of genes detected in the intestine among different challenges, a heat map of RPKM-normalized transcript counts for the top 60 DEGs in each pairwise comparison was generated through hierarchical clustering analysis ( Supplementary Fig. S4). Interestingly, these genes were generally down-regulated after challenges. For example, after LPS injection, we found that 57 of 60 genes were down-regulated. The down-regulated DEG category includes some important elements in immunity and defence response, such as, Complement C1q-like (4-fold), Cathepsin L (232.7-fold), and Heat shock 70 kDa protein 12A (212-fold). Similar results were also observed for Vibrio infection and during fasting (see details in Supplementary Fig. S4). The downregulation of genes related to immunity and defence response may suggest that stressors depress the functions of some immune-related genes in the intestine. The significantly down-regulated genes after all four challenges may be used as biomarkers for detecting stresses in fish.
The DEGs were mapped to different KEGG pathways. The results of mapping to the DEGs are presented in Supplementary Fig. S5 and Supplementary Table S7a. Ortholog assignment and mapping of the DEGs to the biological pathways revealed that 216 (high salinity) to 238 (LPS) pathways were in response to the 4 treatments, respectively (Table 2 and Supplementary Table  S7a). Many pathways previously found acting in neurodegenerative disorder (e.g. Parkinson's disease, Alzheimer's disease, Huntington's disease), regulation of cellular energy metabolism and biosynthesis (e.g. oxidative phosphorylation, cardiac muscle contraction, protein processing in the endoplasmic reticulum, fatty acid metabolism), and signal transduction (Peroxisome proliferator-activated receptor [PPAR] signalling pathway, chemokine signalling pathway, MAPK signalling pathway) responded to the four stressors, suggesting their important functions in stress response. Interestingly, metabolic pathways contain the most abundant DEGs whereby most of the genes were down-regulated after the four challenges. For instance, glyoxylate/hydroxypyruvate reductase and arylamine Nacetyltransferase were down-regulated under LPS and Vibrio injection; 3-hydroxyacyl-CoA dehydrogenase and ubiquinol-cytochrome c reductase cytochrome b subunit were down-regulated after fasting and high-salinity treatments, respectively. The down-regulation of genes in metabolic pathways suggests that the four stressors inhibit metabolic functions of fishes. Ribosome biogenesis is among the most energy-consuming cellular processes, 27 and it is therefore not surprising that the ribosome pathway is tightly controlled upon salinity change (significantly down-regulated by 7-fold), due to the disrupted homeostasis. Mature ribosomes are rapidly degraded by autophagy upon nutrient starvation in Saccharomyces cerevisiae. 27 To our surprise, the number of up-regulated DEGs in the ribosome pathway, e.g. large subunit ribosomal proteins L14e, L18Ae, L19e and L22e, significantly increased by 23-fold, when compared with the number of downregulated DEGs after fasting, implying that selective biosynthesis of ribosomes is functionally important in response to fasting in fish.

Interaction among DEGs and regulated pathways
in response to different stressors Crosstalks between biotic and abiotic stress signalling have been extensively reported in plants, 28,29 yeast and animals. 30 We observed that a lot of DEGs responded to two or more treatments, which may indicate that gene interactions or shared pathways are involved in these stress responses. A Venn diagram in Fig. 2 describes the overlap among genes differentially expressed after different stress challenges. One hundred and twenty-one genes reacted to both Vibrio infection and fasting, and 37 genes were found differentially expressed in all 4 treatments. For example, aquaporins that was involved in osmoregulation 31 and pathogen infection 32 showed responses to all four challenges. The overlapping of a set of DEGs in response to all the four different stressors To analyse the overlaps among differentially regulated pathways involved in the responses, we further compared the differences of regulated DEGs in each pathway under different conditions ( Fig. 3 and Supplementary Table S7b). We found that the DEGs in 75 pathways showed significant up-regulation, and the DEGs in 125 pathways displayed significant downregulation in response to all challenges. Sixty-two pathways were classified into the down-regulated pathway category after all treatments. These pathways were mainly involved in metabolism, digestion and absorption processes, e.g. vitamin digestion and absorption, mineral, fatty acid metabolism, nitrogen, betaalanine, propanoate, pyruvate, fat, glycerolipid, and signalling pathways, e.g. GnRH, neurotrophin, Fc epsilon RI, TGF-beta, phosphatidylinositol and PPAR. Our data suggest that stresses not only inhibit the functions of metabolic pathways, but also repress signalling pathways responsible for reproductive hormone production and growth. Twelve unique pathways were specifically enriched in the up-regulated category after all treatments. These pathways are mainly associated with disorders, signalling transduction involved in innate defences and related processes encoding antimicrobial peptides and proteins, e.g. measles, pertussis, hepatitis C, lysosome, Toll-like receptor (TLR) signalling pathway, notch signalling pathway, NOD-like receptor signalling pathway and cytokine-cytokine receptor interaction. These results suggest that the four stressors induced the expression of genes that related to innate immune defences. Recent studies on mouse revealed that acute and short-term hormonal stress induced an early increase of immune cells in rats 33 and other stresses significantly enhanced the fish immune response. 34 These results supply evidence supporting that stresses enhanced some immune functions for a short time. 35 However, it is not known why the functions of some immune-related genes were enhanced, while other immune-related genes were repressed, which is worthy to further study. PPAR signalling is known to increase in the liver after fasting 36 and is down-regulated in the rat gastrocnemius muscle 37 and in the human skeletal muscle 38 after a 48-h fasting. In our study, down-regulation of gene expression in the PPAR signalling pathway upon LPS, Vibrio infection, conditions of fasting and high-salinity treatments suggests crosstalks for the PPAR signalling pathway among different conditions.
Challenges of fish with stressors can cause specific changes in gene regulation. About 43.67% (645 of 1477; Vibrio infection) to 58.47% (718 of 1228; high salinity) of the DEGs were differentially expressed in a treatment-specific manner, respectively (Fig. 2). These data imply the importance of their roles in response to specific stressors. For example, we observed Annexin A2 (isotig00837), Monoamine oxidase gene (contig00 009), Butyrophilin subfamily 1 member A1 (isotig

00448) and
Alkaline phosphatase (isotig00361) are unique in response to LPS, infection, fasting and high salinity, respectively. Similar functions have been shown for these genes in previous studies. Annexin participates in the establishment of inflammation and the immune response. 13,39 Monoamine oxidase expression in rats reduced after infection with Pneumococcus. 40 Butyrophilin subfamily 1 member A1 is assumed to contribute to milk secretion and lipid droplets, 41 and Alkaline phosphatase is involved in osmotic stress. 42 In addition, 9 (fasting) to 15 (high salinity) regulated pathways were treatment-specifically up-regulated, and 2 (Fasting) to 9 (Vibrio) regulated pathways were treatment-specifically down-regulated (Fig. 3). These unique pathways in response to different stressors generally contained few DEGs (less than four) (Supplementary Table S7b). Little is known about the expression response for most of these DEGs and pathways in fish during stress. The unique kinetic expression of specific genes and unique pathways implies their importance in host stress responses to specific stressors. It would be interesting to characterize their functions in future, such as how stressors are discriminated by host and the mechanisms involved in the stress response.

The effects of stresses on the immune system of the intestine in fish
Although the intestine is not a major immune organ of teleost fish, it plays an importance role in defence against pathogen invasion. 43 It is generally believed that stresses depress the functions of the immune system. 35 Different stresses may have different impacts on immunity. 44 Our intestine transcriptome data set allowed exploring effects of stress on the immune system of the intestines in fish. We found that many DEGs (i.e. 177, 244, 108, and 123 DEGs in response to LPS, Vibrio infection, fasting, and high-salinity challenges, respectively) were classified into the immune system process. These DEGs including interleukin-1, interleukin-8, transforming growth factor, antimicrobial enzymes (e.g. lysozyme, lectin), MHC-related proteins, immunoglobulin superfamily members and pathways such as lysosome and the TLR signalling pathway responding to infection have been discovered (Supplementary  Tables S6 and S7a). These data suggest a dramatic stress-related immune response to stressors in fish intestines. Interestingly, most immune-related DEGs were down-regulated after challenges (see Supplementary Table S6), suggesting that these stresses partially repress the immune function of the intestine.
Fine-tuning of many immune-related DEGs were observed in response to different stressors. Stress induces an increase of a few inflammatory cytokines, such as C -C motif chemokine 25-like significantly increased after Vibrio infection, fasting, and salinity challenges. However, some other cytokines such as tumor necrosis factor ligand superfamily member 10like, granulins-like, tumor necrosis factor receptor superfamily member 14-like were significantly downregulated after the Vibrio infection. The members in the complement system, helping destroy the pathogens and eliminate the infection, also showed differential regulations. For example, complement C3-like is significantly increased after LPS challenge, which is an expected response to acute bacterial infection. However, complement C1q-like protein 2-like showed significant down-regulation after LPS and Vibrio challenges, but substantial up-regulation after fasting and salinity challenges. Our data suggest complicated regulation of cytokines and complement factors in the fish intestines during different stresses. Meprins generate biologically active IL-1b from its precursor pro-IL-1b and play a critical role in the inflammatory response. 45 We found that meprin A subunit beta-like was significantly down-regulated ( 6 -9-fold) after the Vibrio infection, but increased by 2 -3-fold after fasting and salinity challenges, suggesting its important role in different stress responses. Immunoglobulins play critical role in mucosal immunity, acting as the first line of adaptive humoral immune defence at mucosal surfaces in mammals. We found that immunoglobulin lambdalike polypeptide 1-like was down-regulated (2 -3-fold) after LPS injection and Vibrio infection, but significantly increased by 2-fold after fasting, suggesting differential expression level of the adaptive immune system in fish in response to different stressors. Serum amyloid helps the immune system to recognize invasion by bacteria. 46 In the liver of Atlantic salmon, serum amyloid A, complement factor B, and serotransferrin were decreased in expression following starvation, but serum amyloid A and complement factor increased after infection. 47 In our study, serum amyloid P-component-like was highly induced after LPS injection ( 2.6fold), fasting ( 23-fold), and salinity ( 156-fold) challenges, but down-regulated after Vibrio infection (.100-fold). These results suggest that although there is a global coordination of gene expressions under different stresses, there is a fine-tuning of many immune-related DEGs in response to different stressors. Future studies on DEGs related to immune-responses would help to elucidate the complex alterations of the immunological network under stress.

Structure and functional implication of splice
variants of IFABP genes Splice variants can increase the proteome diversity and cellular function. 48 In the RNA-seq data, 7% genes contain splice variants (data not presented). Of which, IFABP (fatty acid-binding protein, intestinal) genes showed significant responses to different challenges. Since these challenges, such as fasting, seriously 456 Transcriptome in Asian Seabass [Vol. 20, affect the nutritional and health status of fish, we speculated that IFABP genes are important to the growth and development of fish. To explore the functions of the IFABP genes, we cloned the genomic DNA regions of the IFABP genes. In comparison of the genomic DNA regions with the transcript sequences and the prediction of ORF for these genes, we identified four splice variants (a1-a4) for IFABP-a and two splice variants (b1-b2) for the IFABP-b gene ( Fig. 4A and Supplementary Fig. S6).
The length for the possible ORF for each variant ranged from 43 amino acids (a4) to 177 amino acids (a2). All of the predicted protein sequences showed sequence conservation with the reported IFABP genes. To find the association of the genes to growth traits, we developed SNP markers from their genomic sequences and conducted QTL mapping by using an Asian seabass F 2 family with 359 individuals and their growth trait data. 49 IFABP-a and IFABP-b were mapped onto the linkage group LG5 and LG14 of the Asian seabass genetic map, respectively. Interestingly, we found that an IFABPa-SNP1245 marker was located near or on a QTL with a proportion of explained phenotypic variance of 10.2% for body weight that were measured at the age of 9 months. 49 IFABP-a-SNP1245 showed high heterozygosity (CG genotype; 77-90%) in the large body weight group (n ¼ 30) and showed low heterozygosity (20-30%) in the small body weight group (n ¼ 30) (Fisher's exact test: P , 0.001). 49 Both analyses suggested that the SNP mutation IFABP-a-SNP1245 was significantly associated with growth. To furthercharacterize the expression pattern of IFABP splice variants, we performed the qPCR gene expression of the splice variants in the intestines of the Asian seabass samples with extreme growth traits (Fig. 4B). We found that all of the four splice variants of the IFABP-a gene and the two splice variants of the IFABP-b gene were significantly highly expressed in the big size group of fishes (P , 0.01), although there were small expression differences among individuals within same groups. The IFABP genes were abundantly expressed in the fish intestines 50,51 and involved in the uptake of dietary fatty acids and their intracellular transport. 50 In mammals, significant associations were found among fatness, the abundance of IFABP and growth. 52,53 The present findings that growth was related to the IFABP expression in fish were consistent with previous researches. 52,53 The higher ability to harvest dietary fatty acids from food as reflected by higher expression of IFABP genes in the intestines of big size fishes may partially explain their fast growth. Our study underlines the importance of splicing variants of IFABP genes in the growth and stress response of the Asian seabass.

Conclusions
This study represents the first RNA-seq-based gene expression analysis in response to several stressors in cultured food fish species, and supplies large transcriptome data for the addressing of biological questions in fish. We detected a broad representation of stressrelated genes in many pathways in response to different challenges. The overlap among DEGs and regulated pathways after different challenges suggest a global coordination in response to different stressors. The differences in DEGs and regulated pathways in the response to distinct types of stress suggest the fine-tuning of gene regulation during different stresses. In addition, we detected dramatic immune responses of the intestine under different challenges. This study is the first step towards achieving a comprehensive view of the molecular mechanisms underlying stress responses. Further detailed analysis of the functions of DEGs and pathways, as well as their interactions will enhance our understanding of the molecular mechanisms involved in the response of fish to stressors.