Transcriptome analysis gene expression in the liver of Coilia nasus during the stress response

Background The estuarine tapertail anchovy (Coilia nasus) is widely distributed in the Yangtze River, the coastal waters of China, Korea, and the Ariake Sound of Japan. It is a commercially important species owing to its nutritional value and delicate flavor. However, Coilia nasus is strongly responsive to stress, this often results in death, which causes huge losses. In this study, we used next-generation sequencing technologies to study changes in gene expression in response to loading stress and the mechanism of death caused by loading stress in Coilia nasus. Results Using next-generation RNA-seq technologies on an Illumina HiSeq 2000 platform, we assembled a de novo transcriptome and tested for differential expression in response to stress. A total of 65,129 unigenes were generated, the mean unigene size and N50 were 607 bp and 813 bp, respectively. Of the assembled unigenes, we identified 2,990 genes that were significantly up-regulated, while 3,416 genes were significantly down-regulated in response to loading stress. Pathway enrichment analysis based on loading stress-responsive unigenes identified significantly stress related pathways. “Metabolism” and “immunity” were the two most frequently represented categories. In the “metabolism” category, “glucose metabolism” and “lipid metabolism” were major subclasses. The transcriptional expression of rate-limiting enzymes in “glucose metabolism” and “lipid metabolism” was detected by RT-qPCR, all were significantly increased after stress. Apoptosis associated proteins tumor necrosis factor alpha (TNF-α), caspase 9, cytochrome c and caspase 3 in the stress group were significantly elevated, moreover, liver injury indicators (Alanine aminotransferase, ALT, and aspartate transaminase, AST) were also significantly elevated, which indicates that loading stress induced liver injury. Conclusion This study provided abundant unigenes that could contribute greatly to the discovery of novel genes in fish. The alterations in predicted gene expression patterns reflected possible responses to stress. Loading stress may induce liver injury through the mitochondrial apoptosis pathway, which was activated by TNF-α. Taken together, our data not only provide information that will aid the identification of novel genes from fish, but also shed new light on the understanding of mechanisms by which physical stressors cause death in fish. Electronic supplementary material The online version of this article (doi:10.1186/1471-2164-15-558) contains supplementary material, which is available to authorized users.


Background
Stress exists in all aspects of aquaculture, daily management such as netting, loading and transporting can lead to a stress response. Stress may result in inhibition of growth, reproductive failure, and reduced resistance to pathogens [1][2][3][4][5]. The estuarine tapertail anchovy (Coilia nasus, junior synonym C. ectenes) is widely distributed in the Yangtze River, the coastal waters of China, Korea, and the Ariake Sound of Japan [6]. It is a commercially important species owing to its nutritional value and delicate flavor. However, Coilia nasus is strongly responsive to stress, this often results in death, which causes huge losses. Currently, little is known of the mechanism by which loading stress causes death in fish. In this study, we used next-generation sequencing technologies to study gene expression changes in response to loading stress and the mechanism of death caused by loading stress in Coilia nasus.
Recently, several species-specific cDNA microarrays have been developed for teleosts, and they are being used increasingly to reveal global gene expression patterns in response to stressor exposure and/or hormone treatment [7][8][9][10]. These studies in fish have revealed that the majority of changes, inspite of the use of different arrays and different types of acute stressor, can be broadly categorized functionally into genes encoding proteins involved in metabolism, immune function and reproduction. Although these studies have identified several genes that were previously not known to be stress responsive, the significance of these observed transcript changes to overall stress adaptation is far from clear.
While the previous studies that used various microarray platforms in fish were limited by large representation of expressed sequence tags (EST), which made gene-specific interpretation of the data difficult, the advent of next-generation sequencing methods, including pyrosequencing, has effectively overcome this limitation [11,12]. Rapid progress in next-generation sequencing technologies has allowed large-scale efficient and economical production of ESTs. Transcriptome sequencing facilitates functional genomic studies, including global gene expression, novel gene discovery, and assembly of full-length genes [13][14][15][16]. To date, available molecular information on Coilia nasus remains sparse, which limits research on the mechanism by which loading stress causes death. The powerful new technologies provide a new opportunity for studies of species without genome reference databases, and non-model organisms.

Experimental animals
Coilia nasus (average weight, 9.6 g) were adapted to the conditions in a 7.0 × 5.0 × 1.0 m 3 aquarium with a water temperature of 24.5 ± 1.0°C, pH 7.2, and dissolved oxygen concentration of 9.2 ± 0.5 mg O 2 /L dechlorinated and aerated water. The fish were fed twice daily, at 7:00 AM and 5:00 PM. At the onset of the experiments, all fish appeared healthy.

Stressing experiment
In March 2012, three 7.0 × 5.0 × 1.0 m 3 ponds were stocked with 120 juvenile Coilia nasus each. The fish were acclimated to the ponds for approximately 14 months before the experiment, at which time the fish were 15 months old. Excess fish were stocked in the ponds in order to ensure subsequent access to the intended numbers. In detail, Coilia nasus prefer to swim towards lighted areas in dark situations. To take advantage of this phototaxis, a sifter was put into the water and lit, following which the fish moved into the sifter. The fish were euthanized immediately with 70 mg/L buffered tricaine methanesulfonate (MS-222). Using this method, five fish were removed from each pond, these 15 fish were the non-stressed controls and were processed immediately (see tissue sampling below). Subsequently, another five fish were netted from each pond and loaded into 75 × 55 × 33 cm3 tanks. After 0.5 h, the fish were euthanized as described above. Compared with the stressed fish, they did not experience the netting and handling used for the control fish. The mean length was 136.98 mm ± 9.26 SEM and mass was 8.86 g ± 1.76 SEM for all fish (n = 30) sampled in this experiment.

Tissue sampling
During the experiments, euthanized fish were submerged immediately in crushed ice to retard degradation of RNA. Blood was collected into ammonium-heparinized capillary tubes after severance of the caudal fin. All fish appeared healthy during dissection and their livers were removed and placed in liquid nitrogen. Plasma was separated by centrifugation. Plasma and liver samples were stored at −80°C until later analysis.
Animal welfare and experimental procedures were carried out in accordance with the Guide for the Care and Use of Laboratory Animals (Ministry of Science and Technology of China, 2006), and were approved by the animal ethics committee of the Chinese Academy of Fishery Sciences.

Cortisol and glucose
Blood glucose measurement was performed by glucoseoxidase-superoxide enzyme endpoint colorimetry, the test kit was obtained from Shanghai Biological Product Research Institute, Ministry of Public Health, China. Cortisol was measured by radioimmunoassay in accordance with the method described by Pickering and Pottinger [17], the test kit was purchased from Beijing Beifang Biotech Research Institute, China. Plasma samples were measured using a Beckman Cx-4 spectrophotometer (Beckman Coulter, Fullerton, CA, USA). Student's t test was used to analyze differences among all treatments (P < 0.05).

RNA sequencing, assembly and annotation
Transcriptome sequencing was carried out on an Illumina HiSeq 2000 platform that generated approximately 100-bp paired-end (PE) raw reads (BGI, Shenzhen, China). After removing adaptor sequences, ambiguous 'N' nucleotides (with the ratio of 'N' greater than 5%) and low quality sequences (with quality score less than 10), the remaining clean reads were assembled using trinity software [18] as described for de novo transcriptome assembly without a reference genome. For homology annotation, non-redundant sequences were subjected to public databases, including NCBI (http://www.ncbi.nlm.nih.gov/) non-redundant protein (Nr) and non-redundant nucleotide (Nt), Swiss-Prot (http://www.ebi.ac.uk/uniprot/), Gene Ontology (GO, http://www.geneontology.org/), Clusters of Orthologous Groups (COG, http://www.ncbi.nlm.nih. gov/COG/) and the Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.jp/kegg/). If the results from different databases conflicted, a priority order of alignments from Nr, Nt, KEGG, Swiss-Prot, GO and COG databases was followed. Comparison with the Nr, Nt and Swiss-Prot databases was carried out using the BlastX algorithm with an E-value cut-off of 0.00001; GO terms at the 2nd level were used to perform GO annotation; COG and KEGG classification were performed using BlastX with an E-value cut-off of 0.00001.

Analysis of differentially expressed genes
To analyze stress-responsive, differentially expressed genes in Coilia nasus, the number of reads for each of the contigs from the two samples was converted to reads per kilo base per million (RPKM) [19]. Following this, the MA-plot-based method with Random Sampling model (MARS) in the DEGseq package was used to calculate the expression abundance of each contig among the analyzed samples. We used an FDR (false discovery rate) to determine the P-value threshold. An FDR < 0.001 was considered to indicate significant expression abundance.
Pathway enrichment analysis identifies significantly enriched metabolic pathways or signal transduction pathways in differentially expressed genes by comparing them with the whole genome background. Bonferroni adjustments [20] were used to estimate levels of significance. After correction for multiple testing, we chose pathways with a Q-value ≤ 0.05 to represent those significantly enriched in differentially expressed genes. The Q-value is defined to be the FDR analog of the P-value. The Q-value of an individual hypothesis test is the minimum FDR at which the test maybe reported as significant.

Gene expression validation
Genes identified in this transcriptome sequencing analysis were validated and quantified by real-time PCR (RT-qPCR). The primers (Additional file 1: Table S1) were designed according to Illumina sequencing data with Primer Premier 5. The prepared total RNA used in RT-PCR analysis was isolated from the same sample as that used for Illumina sequencing. The RT-qPCR was performed on the ABI 7500 real-time PCR system (ABI, USA) using 2× SYBR green real-time PCR mix (Takara, Japan). The PCR amplification was performed in triplicate, using the following cycling parameters: 94°C for 2 min, followed by 40 cycles of 15 s at 94°C, 15 s at 60°C, and 34 s at 72°C. All samples were analyzed in triplicate and the expression of target genes was calculated as relative fold values using the 2 −△△CT method.

Analysis of liver lipid peroxides, tumor necrosis factor alpha and apoptosis-associated proteins
Hepatocyte damage following stress was assessed by measuring alanine aminotransferase (ALT) and aspartate transaminase (AST) activities in plasma, using corresponding detection kits (Jusbio, Shanghai, China) according to the manufacturer's instructions. TNF-α, cytochrome c, caspase-9, and caspase-3 in hepatic tissue were analyzed using an enzyme-linked immunosorbent assay kit (Zhaorui, Shanghai, China), as described by the manufacturer. Lipid peroxides (LPO) were detected by colorimetric methods, using a kit (Jiancheng, Nanjing, China) according to the manufacturer's instructions.

Statistical analysis
Student's t test or ANOVA, where appropriate, were used to identify significant differences among the treatments at p = 0.05. All data in this study were expressed as mean ± standard (S.D.).

Generation and de novo assembly of Coilia nasus transcriptome data
In this study, we performed transcriptome sequencing of two libraries from liver samples in Coilia nasus via an Illumina HiSeq 2000 platform sequencer: 64.8 and 67.8 million reads were obtained from the two libraries. After removing low-quality reads, short reads and reads belonging to mitochondria, a total of 111,053,176 clean reads corresponding to mRNAs were obtained, these reads covered a total of 9,994,785,840 bases (Table 1).
Using the Trinity assembly program, we generated a total of 65,129 unigenes ( Table 2). The length distribution of unigenes larger than 200 bp is shown in Figure 2. The mean unigene size and N50 were 607 bp and 813 bp, respectively. About half of the unigenes (30,582; 47.0%) were ≥ 500 bp and 648 unigenes were > 3,000 bp in length. The largest unigene was 10,911 bp in length (Table 2).
The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis revealed that diverse pathways were represented in our transcriptome dataset: 25,188 unigenes (Additional file 2: Table S2) were assigned to 259 specific pathways. Among them, "metabolism pathways", "regulation of actin cytoskeleton", "pathways in cancer", "focal adhesion", and "endocytosis" were the five most frequently represented pathways (Additional file 2: Table S2). Some important pathways involved in signal transduction were also identified, including  Max length (bp) 10,911 Figure 2 Length distribution of unigenes.

Loading stress-responsive unigenes in Coilia nasus
Unigene expression was estimated by the FPKM method, and differentially expressed genes were identified by referenceto Audic [21]. We found that 2,990 genes were significantly up-regulated, while 3,416 genes were significantly down-regulated in response to loading stress ( Figure 5). The up-regulated and down-regulated genes are listed in Additional file 3: Table S3. The indicated genes included metabolic genes, enzymes, and other immunerelated genes, such as the GTPase gene, threonine-protein kinase and MHC class I heavy chain gene. These genes showed different expression patterns after stress, which implies that they may play an important role in physiological processes associated with stress. Pathway enrichment analysis based on loading stressresponsive unigenes identified significantly stress-related pathways. Among these pathways, "metabolism" and "immunity" were the two most commonly represented categories. In the "metabolism" category, "glucose metabolism" and "lipid metabolism" were major subclasses, and included "glycolysis/gluconeogenesis", "starch and sucrose metabolism", "glycerophospholipid metabolism" and "fat digestion and absorption" pathways. In the "immunity" category, "phagosome", "HTLV-I infection" and "amoebiasis" were the most frequently represented pathways (Table 4).
Lipid peroxides (LPO) are produced in the fatty acid beta-oxidation process, and are cytotoxic. The LPO  content in the stressed group was significantly higher than that in the control group (P < 0.05; Table 5). We measured hepatic TNF-α levels in both control and stressed groups. As shown in Table 5 Plasma aminotransferases were assayed to evaluate hepatotoxicity in Coilia nasusu. The ALT and AST concentrations in the control group were 4.10 ± 0.85 and 41.45 ± 7.79 IU/L, while in the stressed group the concentrations were 8.76 ± 1.08 and 202.12 ± 10.01 IU/L, respectively. The results indicate that significantly increasing plasma levels of ALT and AST developed following stress (Figure 7).

Discussion
The physiological and biochemical response to exposure to stressors has been studied widely in teleosts. One of the evolutionarily conserved stress responses is the rapid elevation in plasma cortisol and plasma glucose levels in response to stressor challenges [22]. In teleosts, the production of cortisol is under the control of the hypothalamus-pituitary-adrenal (HPA) axis [23,24]. Cortisol is released from the adrenal tissue, resulting in up-regulation of the plasma glucose level. Therefore, plasma cortisol and plasma glucose can be used as indicators of the stress response. In our experiment, these parameters were significantly elevated in the stressed fish when compared with the control fish (Figure 1), this suggested that our sample could be used for transcriptome sequencing.
We found that 2,990 genes were significantly upregulated, while 3,416 genes were significantly downregulated in response to loading stress ( Figure 5). On pathway enrichment analysis, "metabolism" and "immunity" were the two most frequently represented categories. In the "metabolism" category, "carbohydrate metabolism" and "lipid metabolism" were major subclasses,    Protein export ko03060 10 3.21E-02 Figure 6 Glucose-6-phosphatase, glucokinase, hormone-sensitive lipase and carnitine actyltransferase I mRNA expression levels in liver in response to stress. The vertical axis shows the relative gene expression levels (means ± SD), bars represent the mean for each replicate ± SD. Replicates or treatments with different letters are significantly different at P < 0.05. and included the "glycolysis/gluconeogenesis", "starch and sucrose metabolism", "glycerophospholipid metabolism" and "fat digestion and absorption" pathways. These carbohydrate metabolic pathways are mainly involved in blood glucose elevation (gluconeogenesis) and glucose utilization (glycolysis). The transcriptional expression of the rate-limiting enzymes (glucose-6-phosphatase and glucokinase) in this process was detected by RT-qPCR, both enzymes were significantly increased following stress ( Figure 6). In addition, plasma glucose levels were significantly elevated after stress ( Figure 1B), which provided further confirmation of the findings. These results are in agreement with studies that have reported higher activities of glycolytic enzymes after exposure to an acute stressor in fish, this may be necessary to cope with the increased energy demand of the liver, including enhanced gluconeogenesis, that is required to re-establish homeostasis [24]. Those two lipid metabolism pathways are involved mainly in fat digestion, absorption and oxidation. Transcript levels for the ratelimiting enzymes (hormone-sensitive lipase and carnitine actyltransferase I) in adipokinetic activity and the fatty acid beta-oxidation process were also significantly increased after stress. Taken together, these results suggest that molecular regulation of enzymes critical for energy substrate mobilization and utilization is a key mechanism involved in acute stress adaptation in fish. However, problems may be associated with these processes, for example, the fatty acid beta-oxidation process may produce lipid peroxides (LPO) [25]. High levels of LPO or plasma glucose may induce an increase in TNF-α [26]. It is now well accepted that trimerization of the respective receptor by TNF leads to the assembly of the death-initiating signaling complex (disc) [27]. The disc consists of the death domain of the receptor, one or several associated proteins, and procaspase-8. Upon assembly of the disc, procaspase-8 is autocatalytically cleaved and activated [27]. The active caspase-8 can either process downstream effector caspases such as caspase-3, −6, and −7 directly [28,29] or activate mitochondria as an amplification mechanism [30]. If mitochondria are involved, caspase-8 cleaves bid, a protein of the Bcl-2 family [31]. The C-terminal portion of the bid molecule inserts into the outer mitochondrial membrane and induces the release of cytochrome c into the cytosol [31,32]. Cytochrome c, together with dATP, apoptosis activating factor-1 (Apaf-1) and procaspase-9, forms the apoptosome, which leads to the activation of caspase-9 and subsequently caspase-3 [33]. Caspase-3 and other effector caspases cleave death substrates, leading to apoptosis [34]. In our study, the levels of TNF-α, caspase 9, cytochrome c and caspase 3 in the liver were significantly elevated in the stressed fish compared with the control fish (Table 5). ALT and AST are important indicators of liver injury, and were also significantly elevated (Figure 7). These lines of evidence suggested that loading stress induces liver cell apoptosis mediated by TNF-α, which causes liver injury.

Conclusions
In the present study, using Illumina sequencing and bioinformatics analysis we analyzed the liver transcriptome of Coilia nasus that were stressed by loading. The main objective of this study was to annotate genes from this transcriptome analysis and identify potential stress

Figure 7
Liver injury induced by stress in Coilia nasus (n = 15). Plasma samples were collected 0.5 h after challenge, and plasma transaminase (ALT and AST) concentrations were determined. Bars represent the mean for each replicate or treatment ± SEM. Replicates or treatments with different letters are significantly different at P < 0.05.