Generation of a reference transcriptome for evaluating rainbow trout responses to various stressors

Water Quality Parameters

Water quality parameters including temperature, DO, pH, CO₂, and ammonia were measured at 0, 3, and 6 hours during the challenges and as summarized in Additional File 1 (Table S1 Water Quality Parameters). The desired parameters for each treatment were obtained within three hours of the beginning of the respective experiments. Water pH values varied between 7.5 and 7.9 for both treatment and control conditions. Concentrations of nitrogenous compounds levels were also within acceptable limits. In the "high temperature" treatment tanks, temperatures were elevated from 14.3°C to an average of 22.6°C within three hours (time 3) and remained elevated until the six hour (time 6) sampling period, at which temperatures averaged 24.8°C. For the "low-temperature" treatment, water temperature was reduced to 1.8°C within three hours (time 3), and was maintained at an average of 2.3°C until time 6. DO levels are one of the most critical water quality parameters in trout culture and exposure to low levels of DO (5.0 - 6.0 mg/l) can result in mortality [54]. However, critical levels are directly correlated with other water quality parameters, like water temperature. DO levels oscillated between 7.16 and 22.40 mg/l across the experimental tanks. Although the average DO concentration at sampling time (time 6) was lower in the re-use water treatment tanks than the control tank (7.8 versus 10.0 mg/l), it remained within acceptable limits for trout production. However, CO₂ levels were significantly higher (67.8 mg/l average) in treatment tanks compared to control (21.7 mg/l). Recommended upper limits for dissolved CO₂ range between 10 and 20 mg/l [55].

Physiological responses

We measured plasma cortisol as an indicator of primary response and glucose, lactate, chloride and lysozyme as indicators of secondary stress responses to detect if stress responses occurred and to characterize differences in responses that are unique to each stressor or shared among stressors. Figures 1, 2 and 3 show the response of plasma glucose, lactate, chloride and cortisol concentrations in addition to lysozyme activity, respectively, compared across all treatments and controls. Values from all control tanks were pooled and presented as a single mean, which was compared to the mean for each stress treatment. Plasma glucose concentrations were significantly lower than control in high temperature and re-use water treatments, and significantly higher than control in the low-temperature treatment. Plasma lactate concentrations increased greater than 2-fold in high-temperature and saline treatments, and also increased after the handling treatment. In contrast, low temperature and re-use water treatments decreased plasma lactate concentrations. High temperature and saline treatments increased, while low temperature and re-use water treatments decreased plasma chloride concentrations compared to control. Lysozyme activity was slightly decreased in the re-use water treatment. Finally, plasma cortisol concentrations increased dramatically, between 5- to 12-fold, compared to the controls, across all treatment groups. The magnitudes and directions of the differences of these parameters between treatments and the control groups reflect unique physiological response portfolios that result from each stressor, such that a sequencing a collection of transcripts isolated from tissues from stressed fish would serve as a qualitative representation of various physiological responses.

Analysis of gene expression

Stress and its resulting increase in cortisol levels have been reported to affect the onset of the apoptosis mechanism [50–52]. The expression of several apoptosis and stress-related candidate genes were analyzed in gill tissue using quantitative real-time PCR (qPCR) to determine if changes in expression of relevant genes occurs during the stress response.

The caspases are proteolytic enzymes which are the core components of the intracellular apoptosis mechanism. Caspases often function in cascades, where an upstream caspase is activated by its interaction with a caspase adaptor [56, 57]. Caspase 8 (casp8) is one of the upstream initiators known to trigger apoptosis, while caspase 3 (casp3) is activated downstream. Figure 4a represents the differential expression (in fold change) of casp3 and casp8 detected with qPCR. Compared to control, the expression of casp3 was up-regulated in the high temperature treatment only. Expression of casp8 increased in the handling, high-temperature, and salinity treatments.

QPCR results for tumor suppressor gene p53 (p53), nuclear protein 1 (nupr1), and heat shock protein 70 b (hsp70b) are reported in Figure 4b. p53[58] has multiple functions, and can be induced by several stressors through different translation mechanisms; its major functions are to regulate apoptosis and growth arrest [59]. The expression of p53 was reduced in the reuse water and salinity treatments and increased in the high temperature treatment (Figure 4b). Nupr1 is also a multifunctional protein and similar to p53, this protein can be induced by different stressors and takes part in apoptosis regulation [60, 61]. Momoda et al.[35] studied gene expression in rainbow trout liver after handling stress, and documented an increased expression after 3 hours of stress and even higher after 21 h. However, in the current study, nupr1 expression was not significantly different from the control for any of the stress treatments. Heat shock proteins (HSP) were initially identified as being induced after stressful thermal conditions; however, it is presently known that this group of proteins also responds to other environmental stressors [37, 62]. Expression of hsp70b was up-regulated in high temperature and salinity treatments, with a dramatic 600-fold increase in expression with the high-temperature stress.

To compensate for osmotic gain of fluid and diffusive loss of ions, freshwater teleosts developed an active ion transport mechanism. Ions like Na⁺, Cl^-, Ca²⁺ are actively transported across the gill epithelium to maintain high levels of ions in the blood. Euryhaline fishes that move from freshwater to salt water transform their gills from ion absorbing mechanisms to an ion-secreting epithelium, resulting in an increase of Na+/K+-ATPase (nakatp) gene expression [63, 64]. A previous study of nakatp expression in rainbow trout gill [64] found that isoforms nakatpα1c and nakatpα3 did not respond to salinity exposure, while isoform nakatpα1a decreased expression within one day after being transferred to salt water and isoform nakatpα1b gradually increased its expression. In this study, expression of the nakatpα1a, nakatp α1b, and nakatpα1c isoforms were significantly reduced by high temperature treatment only. Isoform nakatpα3 showed increased expression in fish exposed to the salinity treatment (Figure 4).

Transcriptome sequencing and assembly

Transcriptome characterization

For the Newbler assembly, only isotigs that are longer than 100 bp were selected for annotation. To avoid redundancy, all isotigs with 90% or higher identities were grouped in one cluster and only a representative from each group was selected for annotation. This step reduced the number of sequences from Newbler assembly from 83,935 to 62,071. Moreover, the 143,560 reads listed on MIRA3 results as debris (singletons), were masked with RepeatMasker using the INRA (National Institute for Agricultural Research, France) Rainbow Trout repeat and cGRASP Salmon masker databases, and 39,210 masked sequences that have more than 300 base pairs of continuous non repetitive segments were also included in the annotation process. After adding the 37,143 MIRA3 contigs; 138,424 sequences were included used for annotations.

Overall, 56,991 sequences were detected as BLASTx hits. The top 10 distributions of BLASTx species are shown in Figure 5, in which Zebrafish (Danio rerio) is the most counted species, followed by Atlantic Salmon (Salmo salar), Pufferfish (Tetraodon nigroviridis), and Rainbow Trout (Oncorhynchus mykiss). In the mapping step, GO annotation data for 50,991 sequences were retrieved from the GO Database. Default settings for Blast2Go annotation was used in the final annotation step and 43,382 sequences were annotated, among them, 29,856 are from Newbler isotigs/contigs, 5,608 are from MIRA3 contigs, and 7,918 are from the masked reads.

Current rainbow trout transcriptome sequence datasets available in NCBI include 287,967 EST sequences in dbEST and the 454-pyrosequencing data submitted to SRA, which include 1,298,911 sequences of a Swanson double-haploid fish and 1,507,274 other rainbow trout cDNA sequences submitted by INRA. In addition, there are 116,463 sequence clusters generated from 391,356 rainbow trout mRNA and EST sequences in NCBI UniGene database and 161,818 rainbow trout transcriptome reference sequences generated by Salem et al.[31]. We compared our sequence reads with all of these datasets, and found that 69.1%, 78.5%, and 56.2% of the reads can be mapped to the rainbow trout sequences in NCBI EST, NCBI SRA for double-haploid fish, and NCBI SRA from INRA, respectively. Moreover, only 67.5% and 55.2% of the reads can be aligned or partially aligned to the rainbow trout cluster sequences in NCBI UniGene database, and the transcriptome reference sequences of Salem et.al., respectively. These comparisons suggest that significant amounts of new data have been generated in our dataset.

SNP Discovery

Four full sib families were used in this experiment, therefore assembly of the transcriptome data was expected to result in the identification of putative SNPs. A total of 57,207 SNPs were identified from the Newbler assembly which represents 30,293 contigs assigned to 21,335 isogroups. The average number of putative SNPs in a contig was 1.9 and in an isogroup 2.7. The maximum number of putative SNPs in a contig was 20, and 15,667 contigs had only one putative SNP. On average, one putative SNP was found every 483 bp. The maximum number of putative SNPs in an isogroup was 392, and 13,254 isogroups had at least two SNPs. Similarly, an additional 808 putative SNPs were identified from the MIRA3 assembly. Overall, the putative SNPs would include 37,358 transitions and 20,657 transversions, similar to the 2:1 ratio observed in some but not all species [65]. To determine the potential for sequence variation to influence gene function we identified the location of putative SNPs in relation to gene sequences based on BLASTx alignments. Overall, 24,479 putative SNPs were determined to fall within exons, 4907 of which were located in the first position of a codon and 4110 the second, both sets expected to impact amino acid sequence and potentially gene function. We also predicted that 15,321 putative SNPs are located in the third position, which as expected is much more than the first and second positions combined. An additional 1322 SNPs were determined to be within exons however discrepancies of SNP positions across multiple alignments resulted in uncertainty of codon position.

It is likely that some of the putative SNPs are actually paralogous sequence variants (PSVs), and therefore are the result of gene duplication events and not sequence variation [66]. In salmonids we expect a high frequency of PSVs due to an evolutionarily recent genome duplication event [67]. However, as all the sequences are transcribed differences may still be relevant to gene function as evolutionarily diverged paralogs sometimes have unique expression profiles [68]. These putative SNPs will be incorporated into a data set previously isolated from a genomic DNA reduced representation library [69] for genetic analyses in rainbow trout (Additional File 3 Tables S3 and S4). Although it would be interesting to evaluate these SNPs for allelic specific expression, our combining multiple tissues from multiple individuals for the creation of a normalized library for this experiment would not produce meaningful quantitative data. However, additional studies of these samples using RNA-Seq will permit such an analysis, yet we must first determine how to distinguish differential expression of allelic variants from that of closely related genes resulting from the genome duplication events.

Fish

Four unrelated families from the NCCCWA even year broodstock population under selective breeding for growth [70] were used for each of the five stressors and control treatments. The fish were tagged and two fish from each family were placed in each of four 193 liter tanks on flow-through water. For each stressor the four tanks included three replicates and one control; a total of 160 fish were used. The experimental fish were acclimated during a two week period at NCCCWA standard culture conditions (13.7°C; CO₂ 25.0 mg/l; DO 15.1 mg/l; salinity 0 ‰). The fish were weighed at the end of the trial, and had average weight of 57.7 g (range: 28.0-99.0 g). Water quality parameters (DO, CO₂, temperature and nitrogenous compounds) were monitored during the acclimation period and at the beginning (time 0), middle (time 3) and end (time 6) of the stress treatments.

Stress challenges

After acclimation, water conditions in each tank were brought to the experimental parameters over the course of three hours (time 3) and then maintained at stress conditions for three hours. The following stress conditions were used: 1) high temperature 25°C; 2) low temperature 2°C; 3) salinity 32‰ (addition of salt, NaCl); 4) Re-use water, presenting lower DO and higher CO₂; and finally 5) handling/crowding stressors, following procedures previously described by Weber and colleagues [71] as adapted from Pottinger et al.[44]. Following each stress treatment, all eight fish from each tank were euthanized by administering a lethal dose of MS-222 and tissues (blood; brain; gill; heart; kidney; liver; spleen and muscle) were collected immediately and frozen in liquid nitrogen.

Plasma metabolite measurements

To identify the physiological responses to stress, plasma metabolite measurements including concentrations of plasma cortisol, glucose, lactate, and chloride in addition to lysozyme activity were measured. Whereas cortisol is a measure of the primary stress response, the others are measures of secondary stress responses. Glucose and lactate serve as measures of metabolic response, chloride as a measure of the osmoregulatory response, and lysozyme activity as a more direct effect of the stressor on the immune system.

Plasma cortisol was measured by radioimmunoassay following procedures described by Redding et al.[72], as modified by Feist et al.[73], but using cortisol antiserum R4866 (provided by Coralie Monroe, University of California-Davis School of Veterinary Medicine, Department of Reproduction), which was validated for use in rainbow trout (Barry et al.[74]). Plasma samples were exposed to heat to denature binding proteins. Glucose was measured by enzymatic coupling with hexokinase and glucose 6-phosphate dehydrogenase (Infinity glucose [HK] liquid stable reagent, Thermo Trace, Melbourne, Australia). Lactate was measured by enzymatic conversion of lactate plus NAD to pyruvate plus NADH (lactate dehydrogenase, Sigma, St Louis MO, #L-2500; NAD Sigma #N-7004). Chloride was measured using a Labconco model 4425000 digital chloridometer (Kansas City MO). Plasma lysozyme activity was measured using a lysozyme turbidity assay as described by Muona and Soivio [75] as modified by Vivas et al.[76]. The bacterium used was Micrococcus lysodeikticus (Sigma, St Louis MO, M3770) and lyophilized hen egg white lysozyme was used as standard (Sigma, L-6876).

Quantitative PCR

Housekeeping genes were used to normalize the transcript expression values [Elongation factor (1f1-α) and β- actin (bact) [77, 78]]. The geNorm applet for Microsoft Excel [72] was used to test the gene expression stability for the candidate genes (M = 0.595 for both 1f1α and bact), and their geometric mean was used to generate a normalization factor for each sample. The normalization factor was used to calculate the expression of each gene relative to the geometric mean expression of the selected reference genes. Expression values were calculated as a fold change, relative to the expression of control group, which was set at 1.0. Fold change data were normalized using log₂ transformation prior to statistical analysis.

Significant effects of stress on gene expression were detected using the general linear models procedure using PC-SAS (Version 9.1). In the event of a significant p-value (P ≤ 0.05), differences between treatments were detected using the Fischer's least significant difference (LSD) procedure. Data on metabolites were analyzed using SigmaStat (Version 3.1; Jandel Scientific, San Rafael, CA). Significant effects of stress on metabolite expression were detected by comparing among control and all treatment groups by Kruskal-Wallis One Way Analysis of Variance on Ranks, followed by Dunn's Method for all pairwise multiple comparisons of means. Treatment means were considered significantly different when P ≤ 0.05. All data are presented as means ± SEM (standard error of the mean).

Transcriptome sequencing

To identify the collection of genes expressed during response to stressors, RNA was extracted from gill, brain, liver, spleen, kidney and muscle from one random sample from each treatment tank plus three control samples (total of 108 RNA samples). All RNA samples were diluted in water and pooled together in a 500 ng/μl concentration. Library preparation and sequencing were performed at the Indiana University's Center for Genomics and Bioinformatics. Methodology was based on Meyer et al.[23], with modifications by K. Mockaitis as described by Schwartz et al.[22]. Briefly, cDNA was synthesized from the RNA pool using primers optimized for GS FLX Titanium sequencing and dsDNA was generated by PCR amplification. Then, cDNA was normalized to reduce sequence coverage of highly expressed transcripts. Accordingly, hybridization and double-stranded nuclease (DSN) digestion were used for the normalization. To optimize sequencing efficiency, libraries were titrated by emulsion PCR/bead enrichment, performed as suggested by Roche 454. Transcriptome sequencing was performed using GS FLX Titanium (Roche 454) pyrosequencing technology.

Transcriptome sequence assembly

The generated rainbow trout cDNA sequence reads were quality trimmed and the adaptor's sequences were removed. Subsequently, GS De Novo Assembler (Newbler) v 2.5p1 was used for sequence assembly. Default parameters for the algorithm of transcriptome assembly were selected, except for specific threshold parameters for generating isotigs. To further ensure the quality of the input reads, the contaminant database included in SeqTrim [79] package for screening of contaminant sequences such as those of Escherichia coli and cloning vectors was used as the optional trimming database of Newbler, in which a read is screened out if it entirely matches with a sequence in the contaminant database. The transcriptome algorithm of Newbler aligns one read completely or partially to a contig. This pipeline also identifies whether a read is repetitive, a singleton or too short (< 50 base pairs) to be used in the computation, or even if it is problematic for assembly. Portions of repetitive reads were included in the assembly results if they aligned to unique contigs. A unique feature of Newbler v2.3 and higher, is the ability of inferring splice-variants in its transcriptome assembly algorithm, in which all contigs are grouped into collections called isogroups based on the contig branching structures and the contig branching structures in each isogroup are further traversed, resulting in a set of isotigs, each representing a possible traversing path of the contigs in the isogroup. Therefore, there is a hierarchical structure among isogroups, isotigs, and contigs. An isogroup, which is analogous to a gene, is composed of a set of isotigs, which can be thought as all possible isoforms (splicing variants) of the gene; and an isotig contains one or more contigs, which can be roughly considered as exons in the gene as demonstrated by Ewen-Campen et al.[80]. In our study, some threshold parameters for generating isotigs were set to the maximum allowed values: the maximum number of isotigs in an isogroup was set to 1000, the maximum number of contigs in an isogroup to 1000, and the maximum number of contigs in one isotig to 200. These parameters enabled the program to generate all possible isotigs for most isogroups. After analysis, several reads were considered singletons by Newbler which included similar sequences, therefore Newbler singleton sequences were re-assembled using MIRA3 (MIRA-3.2.1).

SNP discovery

For SNP identification, high quality sequence reads were first screened for repetitive elements with RepeatMasker [81], using the cGRASP Salmon repeat library [82] and the INRA Rainbow Trout repeat library [83] as reference. Using ssaha2 [84], repeat masked sequences with less than 100 base pairs of consecutive unmasked bases were mapped to the assembly contigs that were longer than 100 base pairs. Then ssaha_pileup pipeline was used to detect sequence variations among the mapped sequences. In order to increase the accuracy of the prediction, the identified putative SNPs were then further processed. First, putative SNPs having total component reads less than five or the component reads for the minor allele less than two were excluded. Then, to address the homopolymer errors commonly found in 454 sequences, only putative SNPs that were at least 50 base pairs apart from other putative SNPs were included. Additionally, flanking sequences were required to be at least 25 base pairs long. Alignment results for the Newbler isotigs and MIRA3 contigs from the BLASTx phase of the BLAST2GO annotation (see below) process were queried to identify query start, query end, and SNP position. These data were used to identify whether or not the SNP occurred inside an exon; if so they were used to calculate if the SNP is in the first, second or third position of the codon. Only data that was in agreement across all BLASTx alignments for each query was considered a positive result.

Sequence Annotation

Functional annotations of Newbler isotigs/contigs, MIRA3 [85, 86] contigs, and unassembled reads were performed using Blast2GO (B2G) [87], a Gene Ontology (GO) based sequence annotation tool. The Blast step was performed using Blastx search against the NCBI non-redundant (nr) protein database, with the expectation value threshold set to 1.0e^-6. Isotigs were first processed with CH-HIT to remove redundant sequences before being used in the annotation step. All isotigs with 90% or higher identities were grouped in one cluster and only a representative from each group was picked for annotation. Functional annotations of Newbler isotigs/contigs and MIRA3 contigs were performed using Blast2GO.

Sequence datasets comparison

A total of 3,091,890 sequences that had been cleaned of vector and E. coli sequences with SeqTrim were compared with the rainbow trout sequences in public available databases using SAHHA2 with the command line parameters: -454 -seeds 5 - score 60 -kmer 13 -skip 4 -diff 0 -output cigar -identity 95. The datasets compared are rainbow trout sequences in dbEST release 080111, sequences in NCBI SRA with the Accessions of SRX007396, SRX041526 - SRX041537, sequence clusters of NCBI UniGene Omy build 32, and the rainbow trout transcriptome reference sequences downloaded from NAGRP Aquaculture Genome Projects (http://www.genome.iastate.edu/aquaculture/salmonids/rainbowtrout/EST_WV.html).