Skip to content

Advertisement

You're viewing the new version of our site. Please leave us feedback.

Learn more

BMC Genomics

Open Access

Transcriptomic response to heat stress among ecologically divergent populations of redband trout

BMC Genomics201516:103

https://doi.org/10.1186/s12864-015-1246-5

Received: 27 January 2014

Accepted: 15 January 2015

Published: 21 February 2015

Abstract

Background

As ectothermic organisms have evolved to differing aquatic climates, the molecular basis of thermal adaptation is a key area of research. In this study, we tested for differential transcriptional response of ecologically divergent populations of redband trout (Oncorhynchus mykiss gairdneri) that have evolved in desert and montane climates. Each pure strain and their F1 cross were reared in a common garden environment and exposed over four weeks to diel water temperatures that were similar to those experienced in desert climates within the species’ range. Gill tissues were collected from the three strains of fish (desert, montane, F1 crosses) at the peak of heat stress and tested for mRNA expression differences across the transcriptome with RNA-seq.

Results

Strong differences in transcriptomic response to heat stress were observed across strains confirming that fish from desert environments have evolved diverse mechanisms to cope with stressful environments. As expected, a large number of total transcripts (12,814) were differentially expressed in the study (FDR ≤ 0.05) with 2310 transcripts in common for all three strains, but the desert strain had a larger number of unique differentially expressed transcripts (2875) than the montane (1982) or the F1 (2355) strain. Strongly differentiated genes (>4 fold change and FDR ≤ 0.05) were particularly abundant in the desert strain (824 unique contigs) relative to the other two strains (montane = 58; F1 = 192).

Conclusions

This study demonstrated patterns of acclimation (i.e., phenotypic plasticity) within strains and evolutionary adaptation among strains in numerous genes throughout the transcriptome. Key stress response genes such as molecular chaperones (i.e., heat shock proteins) had adaptive patterns of gene expression among strains, but also a much higher number of metabolic and cellular process genes were differentially expressed in the desert strain demonstrating these biological pathways are critical for thermal adaptation to warm aquatic climates. The results of this study further elucidate the molecular basis for thermal adaptation in aquatic ecosystems and extend the potential for identifying genes that may be critical for adaptation to changing climates.

Keywords

Climate changeOncorhynchus mykiss gairdneriRedband troutRNA-seqThermal adaptationTranscriptome

Background

Thermal adaptation is a widespread phenomenon in organisms that are exposed to variable and extreme environments. While some organisms may alter their distribution or behavior to avoid stressors and others may acclimate through physiological plasticity [1,2], many species evolve adaptive responses to local conditions over generations through natural selection [3-5]. Evolutionary adaptation to local environments has been demonstrated across a wide variety of taxa [6], and is expected to play a critical role for species with limited dispersal capabilities. However, few studies have identified the underlying molecular mechanisms that have led to conspecific adaptation to thermal conditions.

Molecular techniques such as RNA-seq [7] provide the opportunity to investigate transcriptional response to thermal stress and further identify mechanisms for thermal adaptation. Patterns of gene expression under heat stress are important to determining evolutionary adaptation among conspecific populations that occupy various environments. Multiple genes have been shown to be involved in heat tolerance across many species, including highly conserved heat shock proteins (hsps) that are upregulated under stressful conditions such as exposure to heat [8,9]. An adaptive heat shock response has additionally been shown to occur among conspecific populations that occupy variable environments [3,10]. However, many genes are known to have a role in regulating the effects of temperature and are likely to be involved in thermal adaptation [11,12]. Thus, RNA-seq provides the opportunity to investigate differential expression across the transcriptome and identify biological pathways involved in evolutionary response to thermal stress.

Redband trout (Oncorhynchus mykiss gairdneri) occupy highly variable environments including both montane and desert streams and have been shown to be locally adapted to these different environments [13]. Previous research has demonstrated an adaptive heat shock response among populations from different climates but also suggests that additional mechanisms are involved with thermal adaptation [14]. This species appears to have evolved a finely tuned heat shock response that likely requires additional genes to balance the short term (immediate cellular damage) and long term (fitness) costs associated with thermal stress. Given that oxygen delivery is limiting for fish under climate-related stressors [15], genes involved in oxygen transport are expected to play a significant role. Additionally, we expect that metabolic and immune pathways could be involved given the energy demands and potential for disease under chronic environmental stress [16,17]. Therefore, this study tests for molecular response to heat stress across the transcriptome of ecologically divergent populations of redband trout that have evolved under local climate regimes.

In this study, desert and montane strains of redband trout and their F1 crosses were tested for differential gene expression under heat stress in a common garden experiment. We tested for both acute and chronic stress response by quantifying gene expression in fish that experienced diel water temperatures similar to desert environments that peaked near their thermal maxima (~28.5 C) over the course of four weeks. Tissues were collected from each strain of redband trout at multiple time points (Day 1, 3, 7, and 28) to test for both acclimatization within each strain and evolutionary adaptation between strains. Results are quantified to confirm the role of hsp genes in these fish, but also identify additional genes and biological pathways that are regulated to balance the costs of stress response in populations that have evolved to desert environments.

Results

Sequence alignment to reference transcriptome

A total of 16 lanes for eight pooled libraries provided 2.30 billion quality filtered reads over all 72 samples. Read numbers ranged from 22.40-72.40 M per sample, with an overall mean of 31.96 M (Table 1). Read numbers were well balanced between treatment and control groups with 1.10 billion and 1.21 billion reads, respectively. Trimmed 60 bp reads were aligned with a minimum criterion of 95% identity to the reference transcriptome at an average of 24.1% and 7.53 M reads per sample. Mean percent alignment and mean number of aligned reads for each set of biological replicates ranged from 11.1-40.9% and 3.32-12.49 M, respectively (Additional file 1: Table S1). Reads aligned to a total of 25,128 unique contigs from the reference transcriptome. The mean percentage of aligned reads was not high (24.1%), but this experience is common for non-model species [18]. Trade-offs exist between aligning RNA-seq data to an existing reference transcriptome versus de-novo assembly [7,19,20]. We chose a conservative approach by aligning to a published reference transcriptome [21] that was assembled from multiple types of tissues under stress response which should be more representative than de-novo alignment with a single tissue type.
Table 1

Summary data for redband trout samples including strain (LJ = Little Jacks Cr., K = Keithley Cr., LJxK = F1 crosses), temperature treatment (28°C treatment or 15°C control), sample day, sequencing reads (M = millions), and reference alignment statistics (transcriptome is abbreviated in column heading as Transc.)

      

Mean

 

Mean

    

Quality

Transcr.

Match (M)/

Transc.

Match (%)/

Sample*

Strain

Temp

Day

Reads (M)

Match (M)

Strain

Match (%)

Strain

1

LJ

Trt-28C

1

28.1

4.6

7.0

16.2

23.1

2

LJ

Trt-28C

1

31.3

8.2

 

26.2

 

3

LJ

Trt-28C

1

30.1

8.1

 

26.8

 

4

LJ x K

Trt-28C

1

30.4

12.1

12.5

39.8

40.9

5

LJ x K

Trt-28C

1

27.6

11.6

 

42.1

 

6

LJ x K

Trt-28C

1

33.6

13.7

 

40.9

 

7

K

Trt-28C

1

32.4

13.1

12.1

40.4

40.3

8

K

Trt-28C

1

31.4

13.4

 

42.5

 

9

K

Trt-28C

1

25.9

9.9

 

38.0

 

10

LJ

Trt-28C

3

40.6

14.1

10.8

34.7

30.2

11

LJ

Trt-28C

3

34.8

9.8

 

28.3

 

12

LJ

Trt-28C

3

30.5

8.4

 

27.7

 

13

LJ x K

Trt-28C

3

31.3

9.1

9.0

29.0

26.1

14

LJ x K

Trt-28C

3

30.2

6.9

 

23.0

 

15

LJ x K

Trt-28C

3

41.7

10.9

 

26.3

 

16

K

Trt-28C

3

27.3

9.6

7.0

35.2

24.4

17

K

Trt-28C

3

28.3

7.6

 

27.0

 

18

K

Trt-28C

3

34.3

3.8

 

11.0

 

19

LJ

Trt-28C

7

29.3

12.4

9.1

42.4

29.2

20

LJ

Trt-28C

7

31.5

9.4

 

29.9

 

21

LJ

Trt-28C

7

35.1

5.4

 

15.5

 

22

LJ x K

Trt-28C

7

28.5

10.8

10.5

37.9

38.2

23

LJ x K

Trt-28C

7

26.4

9.9

 

37.6

 

24

LJ x K

Trt-28C

7

28.0

10.9

 

39.1

 

25

K

Trt-28C

7

37.4

6.8

5.2

18.2

16.6

26

K

Trt-28C

7

28.7

5.0

 

17.3

 

27

K

Trt-28C

7

27.8

4.0

 

14.2

 

28

LJ

Trt-28C

28

22.8

6.4

5.6

28.2

20.6

29

LJ

Trt-28C

28

31.2

7.5

 

24.0

 

30

LJ

Trt-28C

28

29.6

2.8

 

9.6

 

31

LJ x K

Trt-28C

28

43.0

0.5

3.7

1.1

14.9

32

LJ x K

Trt-28C

28

24.1

4.1

 

16.9

 

33

LJ x K

Trt-28C

28

24.3

6.4

 

26.5

 

34

K

Trt-28C

28

26.4

6.3

6.0

23.7

23.0

35

K

Trt-28C

28

25.4

7.3

 

28.7

 

36

K

Trt-28C

28

26.8

4.5

 

16.6

 

37

LJ

Con-15C

1

23.8

7.4

6.9

31.2

27.0

38

LJ

Con-15C

1

24.7

6.7

 

27.3

 

39

LJ

Con-15C

1

29.4

6.6

 

22.5

 

40

LJ x K

Con-15C

1

25.0

6.5

10.2

25.9

25.0

41

LJ x K

Con-15C

1

65.0

15.3

 

23.5

 

42

LJ x K

Con-15C

1

34.2

8.7

 

25.6

 

43

K

Con-15C

1

22.4

5.3

6.9

23.7

26.0

44

K

Con-15C

1

27.1

8.3

 

30.6

 

45

K

Con-15C

1

29.6

7.0

 

23.7

 

46

LJ

Con-15C

3

24.2

8.1

5.4

33.5

19.7

47

LJ

Con-15C

3

34.4

3.7

 

10.9

 

48

LJ

Con-15C

3

30.0

4.4

 

14.8

 

49

LJ x K

Con-15C

3

46.3

0.7

3.3

1.6

11.1

50

LJ x K

Con-15C

3

27.0

3.6

 

13.4

 

51

LJ x K

Con-15C

3

30.6

5.6

 

18.4

 

52

K

Con-15C

3

26.7

8.3

5.8

31.0

21.2

53

K

Con-15C

3

28.0

4.3

 

15.4

 

54

K

Con-15C

3

28.5

4.9

 

17.1

 

55

LJ

Con-15C

7

24.3

7.1

6.0

29.4

24.0

56

LJ

Con-15C

7

26.2

7.0

 

26.9

 

57

LJ

Con-15C

7

25.0

4.0

 

15.9

 

58

LJ x K

Con-15C

7

33.2

6.7

4.3

20.2

13.8

59

LJ x K

Con-15C

7

32.6

2.9

 

8.9

 

60

LJ x K

Con-15C

7

27.3

3.4

 

12.4

 

61

K

Con-15C

7

25.6

5.0

10.8

19.4

23.4

62

K

Con-15C

7

27.8

5.9

 

21.0

 

63

K

Con-15C

7

72.7

21.7

 

29.9

 

64

LJ

Con-15C

28

35.0

12.5

9.5

35.6

25.9

65

LJ

Con-15C

28

41.1

7.1

 

17.4

 

66

LJ

Con-15C

28

36.1

8.9

 

24.6

 

67

LJ x K

Con-15C

28

31.5

5.9

8.3

18.6

21.2

68

LJ x K

Con-15C

28

54.7

7.1

 

12.9

 

69

LJ x K

Con-15C

28

37.3

11.9

 

32.1

 

70

K

Con-15C

28

44.4

0.3

4.9

0.7

13.4

71

K

Con-15C

28

41.5

6.9

 

16.5

 

72

K

Con-15C

28

32.3

7.4

 

22.9

 

mean

--

--

--

31.96

7.53

7.5

24.1

24.1

*Each sample includes 3 pooled RNA samples from the same rearing tank.

Principal component analyses of overall gene expression data clearly showed that samples clustered by treatment or control condition and also that distinct clusters were present for the desert and montane strains, but the cluster for the F1 strain overlapped with the montane strain (Figure 1).
Figure 1

Principal components analysis of overall transcriptome expression. Results for 18 samples of redband trout collected either at first exposure to heat stress up to 28°C (darker shades) or on the same day from 15°C control temperature (lighter shades). Samples are color coded by their environment: desert strain (red; LJ = Little Jacks Cr.), F1 crosses (gray; LJxK), montane strain (blue; K = Keithley Cr.).

Differential expression

Differences in gene expression were highly significant within each strain (Figure 2a), with 7,051 significant genes for the desert strain (4,238 upregulated, 2,813 downregulated), 6,906 for F1 crosses (3,375 upregulated, 3,531 downregulated), and 6,774 for montane (3,499 upregulated, 3,275 downregulated). As expected, a large number of total transcripts (12,814) were differentially expressed in the study (FDR ≤ 0.05) with 2,310 transcripts occurring in common among all three strains, but the desert strain had a larger number of unique differentially expressed transcripts (2,875) than the montane (1,982) or the F1 (2,355) strain (Figure 2a). Strongly differentiated genes (>4 fold change and FDR ≤ 0.05; Figure 2b) were highly abundant in the desert strain (824 unique transcripts; Additional file 1: Table S1) relative to the other two strains (montane = 58; Additional file 2: Table S2; F1 = 192; Additional file 3: Table S3), particularly in upregulated transcripts (Figure 3a-c).
Figure 2

Venn diagrams of differentially expressed genes. Results of each strain of redband trout for a) all significant transcripts (FDR ≤ 0.05); and b) strongly differentiated transcripts (>4 fold change and FDR ≤ 0.05). Circles are color coded to represent fish by their environment: desert strain (red; LJ = Little Jacks Cr.), F1 crosses (gray; LJxK), montane strain (blue; K = Keithley Cr.). For a) there were 12,314 genes that were not statistically significant at either level and are listed outside of the circles on the Venn diagram.

Figure 3

Adaptation: differential expression for each strain of redband trout across all time periods of heat stress versus fish held at control temperatures. Results for a) desert strain from Little Jacks Cr. b) F1 crosses; and c) montane strain from Keithley Cr. Genes that are significantly differentiated (FDR ≤ 0.05) are in red and those that are not significant are in black. On a log2 scale, the green lines show genes that are ≥ 2 fold, and the blue lines designate genes that are ≥ 4 fold. The x-axis is the mean expression of each gene in counts per million reads (CPM) on a log2 scale.

Differences in gene expression were also observed over time as the number of significant genes consistently decreased with more days of exposure to temperature stress (Figure 4). The number of significant genes was 7,833 at Day 1 (4,058 upregulated, 3,775 downregulated), 6,408 at Day 3 (3,344 upregulated, 3,064 downregulated; 18.2% decrease from Day 1), 3,624 at Day 7 (1,958 upregulated, 1,666 downregulated; 53.7% decrease from Day 1), and 1,269 at Day 28 (719 upregulated, 550 downregulated; 83.8% decrease from Day 1). This trend was consistent with the expectation that the stress response would become reduced with chronic exposure to heat stress.
Figure 4

Acclimation: differential expression at each time period of heat stress versus fish held at control temperatures across all strains. Results for a) Day 1 of heat stress treatment; b) Day 3 of heat stress treatment; c) Day 7 of heat stress treatment; and d) Day 28 of heat stress treatment. The number of significant genes was 7,833 at Day 1 (4,058 upregulated, 3,775 downregulated), 6,408 at Day 3 (3,344 upregulated, 3,064 downregulated; 18.2% decrease from Day 1), 3,624 at Day 7 (1,958 upregulated, 1,666 downregulated; 53.7% decrease from Day 1), and 1,269 at Day 28 (719 upregulated, 550 downregulated; 83.8% decrease from Day 1). Genes that are significantly differentiated (FDR ≤ 0.05) are in red and those that are not significant are in black. On a log2 scale, the green lines show genes that are ≥ 2 fold, and the blue lines designate genes that are ≥ 4 fold. The x-axis is the mean expression of each gene in counts per million reads (CPM) on a log2 scale.

Gene ontology and enrichment with Blast2GO revealed that strongly differentiated genes in each strain (>4 fold change and FDR ≤ 0.05) included several categories for each of biological processes, molecular function, and cellular components (Figure 5a-c; Additional files 4, 5 and 6: Tables S4-6). Within biological process, there were a total of 18 pathway categories at level 2 gene ontology, but nearly 70% of the genes were included in five categories: cellular process (mean = 17.2%), metabolic process (mean = 16.7%), response to stimulus (mean = 12.2%), single-organism process (mean = 12.4%), and biological regulation (mean = 10.3%; Figure 5a). Within molecular function, there were a total of 11 pathway categories at level 2 gene ontology with over 75% of the genes in two categories: binding (mean = 52.3%), and catalytic activity (mean = 24.3%; Figure 5b). Within cellular components, there were a total of 10 pathway categories at level 2 gene ontology, with over 80% of the genes included in four categories: cell (mean = 31.3%), organelle (mean = 24.8%), membrane (mean = 16.1%), and macromolecular complex (mean = 10.5%; Figure 5c). Finally, many of the genes that were found to be differentially expressed relative to control fish in all strains and time points were stress response genes such as various hsp transcripts.
Figure 5

Gene ontology (GO) annotation for transcripts that were strongly differentiated expressed in each strain (>4 fold change and FDR ≤ 0.05). Results shown for level 2 categories for a) Biological process (“org.” = organization); b) Molecular function (“TFA” = transcription factor activity); c) Cellular component. Bars are color coded to represent fish by their environment: desert strain (red; LJ = Little Jacks Cr.), F1 crosses (gray; LJxK), montane strain (blue; K = Keithley Cr.).

Patterns of gene expression for each strain over time were compared with results from qPCR assays for heat shock genes and were highly consistent with either RNA-seq or qPCR data (Additional file 7: Figure S1). Specifically, expression patterns showed that heat shock genes were significantly lower for the desert strain at Day 1 for all hsp genes and all strains had decreased gene expression from Day 3 through the remainder of the experiment as shown previously [14].

Discussion

Since thermal stress has broad biological effects on organisms, transcriptional response is expected to be highly diverse across several genes in ectothermic species such as redband trout. This study confirms that numerous genes are differentially expressed in redband trout under heat stress, and several pathways are involved. However, there were key pathways that contained a large proportion of differentially expressed transcripts including response to stimulus, metabolic processes, cellular processes, molecular binding function, and cell membrane function. These pathways correspond well with previous studies that demonstrate these as critical physiological components involved with of aquatic ectotherms exposed to elevated water temperatures [17,22]. In particular, several physiological studies have linked thermal tolerance with aerobic scope and emphasize the role of metabolic processes in thermal adaptation (e.g., [22-24]). The larger number of differentially expressed genes in the desert strain versus the other two strains suggests that a complex combination of genes has evolved for redband trout in their desert environment. It is also possible that a few genes with large pleiotropic effects could be responsible for the pattern observed.

Evidence for acclimation to heat stress was extensive as the number of differentially expressed transcripts decreased by 83.7% from Day 1 to Day 28. Results from this study elaborate on previous findings in redband trout that stress response genes are highly upregulated when exposed to heat stress [25,14]. Multiple heat shock genes (e.g., hsp70, hsp90, hsp47) were differentially expressed in all strains and time periods. However, an acclimation effect was evident as expression levels decreased over time in all strains. This is consistent with theories of acclimation to heat stress where organisms are able to moderate their heat shock response over time, as opposed to initial exposure where immediate survival is a priority [26]. However, the current results demonstrate that acclimation within strains occurs throughout much of the transcriptome and is not limited to heat shock genes.

More importantly, this study demonstrates that adaptive patterns of expression have evolved in ecologically divergent populations of this species. Results from Narum et al. [14] specifically highlight the adaptive response of heat shock genes in redband trout, with lower hsp gene expression observed in desert versus montane strains. Results in heat shock genes from the current RNA-seq data corroborate the previous qPCR results and emphasize that warm adapted natural populations are likely to have evolved a specialized heat shock response that reduces physiological costs of hsp production. This result is consistent with the adaptive heat shock response observed in natural populations of other organisms such as killifish (Fundulus heteroclitus; [27]) and Drosophila buzzatii [10]. This remains an important finding of this study and provides clarification regarding evolutionary adaptation of hsp gene expression in heat tolerant populations. However, many recent studies indicate that complex mechanisms are involved in thermal adaptation of aquatic ectotherms beyond heat shock response (e.g., [22,28,29]. Indeed, this study of the transcriptome revealed adaptive patterns in metabolic and cellular process genes that suggest desert fish are more efficient at supporting these pathways than montane fish under heat stress, and chronic exposure may cause failure of these genes to be expressed in montane and F1 crosses and suggests that some critical physiological functions become limited in these strains over time. Previous studies suggest that metabolic pathways may be particularly important since metabolic energy stores are positively correlated with physiological function and swimming behavior in thermally adapted redband trout [16,17]. The general pattern observed across the transcriptome for F1 crosses indicates intermediate gene expression consistent with additive variation, but there may be a maternal or dominant effect at some genes since more differentially expressed genes were shared between the F1s and the maternal montane strain than the desert strain.

A variety of immediate and long-term anthropogenic disturbances such as habitat destruction and climate change have negative impacts on freshwater fish [30,31] and the need to understand mechanisms for thermal adaptation in these organisms is critical. Many fishes have already been extirpated from large portions of their historical range (e.g., [32]) and the effects of climate change are expected to further alter species’ range, phenology, and persistence [32-35]. Genomic and physiological mechanisms for thermal adaptation can be important tools for conservation measures to enable long-term viability of wild populations [36-38]. Specifically, this study helps to further identify genomic tools such as genetic screening with candidate genes that may be integrated with measurements of cardiac function [22] in order to screen broadly across species’ range to predict the potential for adaptation under scenarios of climate change [39].

Conclusions

This study demonstrates that redband trout from a desert climate have a much larger number of strongly differentially expressed genes than montane and F1 strains in response to heat stress, suggesting that a combination of genes has evolved for redband trout to adapt in their desert environment. Recent studies of physiological adaptation in aquatic ectotherms indicate that intraspecific thermal tolerance is set by limitations in aerobic performance, specifically the upper limit of heart rate to deliver oxygen to tissues (e.g., [15,22,23]). This is due to temperature dependent oxygen limitation in aquatic environments, a theory that has been well supported in many organisms [40]. In order to support this increase in cardiac performance, redband trout would need to differentially regulate genes from multiple pathways including those observed in this study (e.g., metabolic pathways). However, further studies that specifically link individual gene expression [41] with physiological functions such as aerobic scope and heart rate are needed to further elucidate the specific mechanisms involved with thermal adaptation in this species. Development of a reference transcriptome and genome that are more specific to this subspecies of O. mykiss would also provide better annotation of genomic architecture of various traits such as thermal adaptation.

Methods

Redband trout populations and thermal stress experiments

To investigate thermal acclimation and adaptation of redband trout (O. mykiss gairdneri) from desert and montane populations, fry of approximately four months of age from each environment and their F1 crosses were exposed to diel temperature cycles (peaking at 28°C) over a 4-week period in a controlled setting. Gill tissues were collected from euthanized individuals on day 1, 3, 7 and 28 to quantify mRNA expression across the transcriptome. Gill tissue was used for this study since oxygen transport and osmoregulation are expected to play important functional roles in thermal adaptation.

Gametes and fry were collected from two ecologically divergent populations: one from a desert climate stream Little Jacks Cr. (LJ), and one from a montane climate stream Keithley Cr. (K), both located in Idaho, USA. In order to create F1 crosses, gametes from each strain were cross fertilized and reared in laboratory incubators. F1s were included to investigate additive genetic variation associated with response to heat stress. The two sites were chosen for study based on previous tests of redband trout from six desert and six montane streams that demonstrated that Little Jacks Cr. fish were adapted to a desert climate and Keithley Cr. was a typical montane stream population [13]. Gametes were fertilized to produce half-sibling progeny representing three distinct strains: one of pure desert strain (LJ), one of pure montane strain (K), and the F1 crosses (LJ males × K females). Fry were reared in constant 15°C spring water until they reached an average weight of 2 g, then each strain was divided into treatment and control groups. Three replicate tanks were used for all treatment and control groups for each strain (3 strains × 2 treatments × 3 replicates equals a total of 18 tanks) with an average of 45 fish per tank. Fish were fed a diet of Soft Moist pellets (Rangen Inc.) to satiation twice per day, and photoperiod was fixed at 14 h light and 10 h darkness. Fish in recirculating treatment tanks experienced diel temperature cycles over 6 weeks that reached a maximum of 28.5°C in the afternoon and a minimum of 17.0°C at night (mean temperature gradient of ~1.5°C per hour; Additional file 7: Figure S1, Supporting information), while fish in control tanks were held at a constant temperature of 15°C spring water. All experimental protocols were approved in advance by the University of Idaho’s Institutional Animal Care and Use Committee (Protocol #201025).

RNA-seq library prep and Illumina sequencing

Total RNA was isolated from approximately 5 mg of gill tissue from individual fish using Qiagen RNeasy kits. RNA was normalized to 100 ng/μL and equal volumes of RNA from three fish from each tank were pooled for a total of 72 libraries (18 tanks × 4 time periods each; Table 1). The Ribo-Zero™ Magnetic Gold Kit (Epicentre) was used to deplete the samples of ribosomal RNA (rRNA) which constitutes a large proportion of the total RNA. Ribosomal RNA depleted samples were purified using the ethanol precipitation method suggested in the Ribo-Zero Kit protocol and resuspended in 20 μL of RNAse free water. RNA-seq libraries were prepared using 4.75 μL of template material using the strand-specific library preparation kit ScriptSeq™ v2 RNA-Seq Kit (Epicentre, Madison WI USA). The tagged cDNA constructs were purified using Qiagen MinElute columns and sample specific index sequences were added during PCR [95°C – 1 m; (95°C – 30s; 55°C – 30s; 68°C – 3 m) × 15; 68°C – 7 m; 4°C – hold] using ScriptSeq (Epicentre) index PCR primers. Each amplified library was then purified using Agencourt AMPure XP magnetic beads and eluted with 20uL of nuclease free water. Dilutions (1:2000) of the purified and indexed libraries were then quantified by qPCR using a ABI-7900 instrument (Life Technologies), Power-Sybr master mix (Life Technologies), standard Illumina (P5, P7) primers and an Illumina PhiX library as a standard. The indexed libraries were normalized to 10 nM concentration in Tris-EDTA (pH 8.0) buffer with 0.1% Tween-20 and combined for sequencing (8 pooled libraries with nine samples, each with ScriptSeq index sequences 1–9). Each of the pooled libraries was sequenced in two lanes of a single read 100 bp flow cell on an Illumina HiSeq 1500 instrument for a total of 16 lanes. Each lane of data was demultiplexed by index sequence and reads were combined from both lanes for each sample. The average number of reads per sample after quality filtering was 31.96 M and ranged from 22.40 – 72.70 M. Raw data was submitted to NCBI’s short read archive (SRA; entry GSE53907).

Sequence alignment to reference transcriptome

Raw sequencing data was aligned to a reference transcriptome for rainbow trout designed for stress response [21] using the program Bowtie [42]. Parameters for Bowtie were set to exclude the first 10 bases and last 30 bases of sequencing data leaving 60 bases of high quality sequence for alignment. Both the forward (+) and reverse complement (−) of each reference transcript were considered for alignment since the reference transcriptome was assembled from non-directional sequence data. The best single match to the reference transcriptome was returned provided it had no more than 3 mismatches across the 60 bases (95% identical).

Output data from Bowtie was condensed by counting matches to each of the reference transcript contigs for both + and - orientation. Since the library preparation process generates directional constructs (strand-specific library preparation), we expected legitimate alignments to match only one of the two orientations of each reference transcript. Indeed, we found that our library reads predominantly matched only 1 of the 2 strand orientations (>90% of reads). Therefore, read alignments to the minor strand were considered mis-assigned and excluded from the data set. Finally, since the reference transcriptome included several sequence variants (contigs) for many putative genes, combined read counts from contigs within each gene were utilized for differential expression analyses.

Differential expression analyses

Tests for differential expression were completed using the edgeR Bioconductor package [43]. Differential expression analyses were done to test for two processes involved in heat stress response: 1) acclimatization to chronic heat stress over time, and 2) evolutionary adaptation of the specific strains. To test for acclimatization, differential expression was tested separately for each time period (Day 1, 3, 7, 28) with all strains in the model. To test for evolutionary adaptation, differential expression was tested separately for each strain (LJ, K, LJxK) with all time periods in the model. Each model included the additional factor of condition (treatment or control) and three biological replicates. Genes that were not expressed in either condition were removed; specifically only genes with at least two counts per million reads in at least nine samples were kept for further analyses. Gene counts were normalized with the trimmed mean of ‘M’ values (TMM) method in edgeR as this has been shown to be one of the most reliable methods for this purpose in RNA-seq studies [44]. Sample metadata and normalized gene expression data was submitted to NCBI’s Gene Expression Omnibus (GEO; entry GSE53907). As suggested by McCarthy et al. [45], genewise dispersion and a general linear model (GLM) were used for tests of differential expression. Genewise dispersion estimates deprioritize genes with inconsistent results and allow the main analysis to focus on changes that are consistent between biological replicates. The GLM accounts for the multifactor design of this study. A false discovery rate (FDR at 0.05; [46]) was applied to account for multiple tests of differentially expressed genes. In order to validate patterns of gene expression in the RNA-seq data, results were compared to quantitative PCR (qPCR) data for multiple heat shock genes with Sybr Green on an ABI 3730 instrument as detailed in Narum et al. [14]. Significantly regulated contigs in each strain were annotated in Blast2GO with a blastx minimum e-value set to 1.0E-06 [47].

Availability of supporting data

Raw data was submitted to NCBI’s short read archive (SRA; entry GSE53907). Sample metadata and normalized gene expression data was submitted to NCBI’s Gene Expression Omnibus (GEO; entry GSE53907).

Declarations

Acknowledgements

We thank Kevin Meyer, Steve Elle, Liz Mamer, Chris Sullivan, Carson Watkins for assistance with field sampling, Lindsay Maier, Mike Casten, Ron Hardy for experimental design and fish rearing, Stephanie Harmon for assistance in the lab, and Mike Miller for sequencing pipeline suggestions. Funding was provided by the Bonneville Power Administration through grant 200900500 to S.R.N. All experimental protocols were approved in advance by the University of Idaho’s Institutional Animal Care and Use Committee (Protocol #201025).

Authors’ Affiliations

(1)
Columbia River Inter-Tribal Fish Commission

References

  1. Angilleta MJ. Thermal Adaptation: a Theoretical and Empirical Synthesis. Oxford: Oxford University Press; 2009.View ArticleGoogle Scholar
  2. Whitehead A, Galvez F, Zhang S, Williams LM, Oleksiak MF. Functional genomics of physiological plasticity and local adaptation in killifish. J Hered. 2011;102:499–511.View ArticlePubMed CentralPubMedGoogle Scholar
  3. Dahlhoff EP, Rank NE. Functional and physiological consequences of genetic variation at phosphoglucose isomerase: Heat shock protein expression is related to enzyme genotype in a montane beetle. Proc Natl Acad Sci. 2000;97:10056–61.View ArticlePubMed CentralPubMedGoogle Scholar
  4. Hoffman AA, Sorensen JG, Loeschcke V. Adaptation of Drosophila to temperature extremes: bringing together quantitative and molecular approaches. J Therm Biol. 2003;28:175–216.View ArticleGoogle Scholar
  5. Kavanagh KD, Haugen TO, Gregersen F, Jernvall J, Vollestad LA. Contemporary temperature-driven divergence in a Nordic freshwater fish under conditions commonly thought to hinder adaptation. BMC Evol Biol. 2010;10:350.View ArticlePubMed CentralPubMedGoogle Scholar
  6. Keller I, Seehausen O. Thermal adaptation and ecological speciation. Mol Ecol. 2012;21:782–99.View ArticlePubMedGoogle Scholar
  7. Wolf JBW. Principles of transcriptome analysis and gene expression quantification: an RNA-seq tutorial. Mol Ecol Res. 2013;13:559–72.View ArticleGoogle Scholar
  8. Morimoto RI, Sarge KD, Abravaya K. Transcriptional regulation of heat shock genes. J Biol Chem. 1992;267:21987–90.PubMedGoogle Scholar
  9. Sørensen JG, Kristensen TN, Loeschcke V. The evolutionary and ecological role of heat shock proteins. Ecol Lett. 2003;6:1025–37.View ArticleGoogle Scholar
  10. Sørensen JG, Dahlgaard J, Loeschcke V. Genetic variation in thermal tolerance among natural populations of Drosophila buzzatii: down regulation of Hsp70 expression and variation in heat stress resistance traits. Funct Ecol. 2001;15:289–96.View ArticleGoogle Scholar
  11. Sørensen JG, Norry FM, Scannapieco AC, Loeschcke V. Altitudinal variation for stress resistance traits and thermal adaptation in adult Drosophila buzzatii from the New World. J Evol Biol. 2005;18:829–37.View ArticlePubMedGoogle Scholar
  12. Kassahn KS, Crozier RH, Ward AC, Stone G, Caley MJ. From transcriptome to biological function: environmental stress in an ectothermic vertebrate, the coral reef fish Pomacentrus moluccensis. BMC Genomics. 2007;8:358.View ArticlePubMed CentralPubMedGoogle Scholar
  13. Narum SR, Campbell NR, Kozfkay CC, Meyer KA. Adaptation of redband trout in desert and montane environments. Mol Ecol. 2010;19:4622–37.View ArticlePubMedGoogle Scholar
  14. Narum SR, Campbell NR, Meyer KA, Miller MR, Hardy RW. Thermal adaptation and acclimation of ectotherms from differing aquatic climates. Mol Ecol. 2013;22:3090–9097.View ArticlePubMedGoogle Scholar
  15. Pörtner HO, Knust R. Climate change affects marine fishes through the oxygen limitation of thermal tolerance. Science. 2007;315:95–7.View ArticlePubMedGoogle Scholar
  16. Gamperl AK, Rodnick KJ, Faust HA, Venn EC, Bennett MT, Crawshaw LI, et al. Metabolism, swimming performance, and tissue biochemistry of high desert redband trout (Oncorhynchus mykiss ssp.): evidence for phenotypic differences in physiological function. Physiol Biochem Zool. 2002;75:413–31.View ArticlePubMedGoogle Scholar
  17. Rodnick KJ, Gamperl AK, Lizars KR, Bennett MT, Rausch RN, Keeley ER. Thermal tolerance and metabolic physiology among redband trout populations in south-eastern Oregon. J Fish Biol. 2004;64:310–35.View ArticleGoogle Scholar
  18. Hornett EA, Wheat CW. Quantitative RNA-seq analysis in non-model species: assessing transcriptome assemblies as a scaffold and the utility of evolutionary divergent genomic reference species. BMC Genomics. 2012;13:361.View ArticlePubMed CentralPubMedGoogle Scholar
  19. De Woody JA, Abts KC, Fahey AL, Ji Y, Kimble SJA, Marra NJ, et al. Of contig and quagmires: next-generation sequencing pitfalls associated with transcriptomic studies. Mol Ecol Res. 2013;13:551–8.View ArticleGoogle Scholar
  20. Singal S. De novo transcriptomic analyses for non-model organisms: an evaluation of m ethods across a multi-species data set. Mol Ecol Res. 2013;13:403–16.View ArticleGoogle Scholar
  21. Sánchez CC, Weber GM, Gao G, Cleveland BM, Yao J, Rexroad CR. Generation of a reference transcriptome for evaluating rainbow trout responses to various stressors. BMC Genomics. 2011;12:626.View ArticlePubMed CentralPubMedGoogle Scholar
  22. Eliason EJ, Clark TD, Hague MJ, Hanson LM, Gallagher ZS, Jeffries KM, et al. Differences in thermal tolerance among sockeye salmon populations. Science. 2011;332:109–12.View ArticlePubMedGoogle Scholar
  23. Pörtner HO. Climate variations and the physiological basis of temperature dependent biogeography: systemic to molecular hierarchy of thermal tolerance in animals. Comp Biochem Physiol. 2002;132A:739–61.View ArticleGoogle Scholar
  24. Farrell AP. Environment, antecedants and climate change: lessons from the study of temperature physiology and river migration of salmonids. J Exp Biol. 2009;212:3771–80.View ArticlePubMedGoogle Scholar
  25. Cassinelli JD, Moffitt CM. Comparison of growth and stress in resident redband trout held in laboratory simulations of montane and desert summer temperature cycles. Trans Am Fish Soc. 2010;139:339–52.View ArticleGoogle Scholar
  26. Somero GN. The physiology of climate change: how potentials for acclimatization and genetic adaptation will determine “winners” and “losers.”. J Exp Biol. 2010;213:912–20.View ArticlePubMedGoogle Scholar
  27. Fangue NA, Nann AF, Hofmeister M, Schulte PM. Intraspecific variation in thermal tolerance and heat shock protein gene expression in common killifish, Fundulus heteroclitus. J Exp Biol. 2006;209:2859–72.View ArticlePubMedGoogle Scholar
  28. Dong Y, Somero GN. Temperature adaptation of cytosolic malate dehydrogenases of limpets (genus Lottia): differences in stability and function due to minor changes in sequence correlate with biogeographic and vertical distributions. J Exp Biol. 2008;212:169–77.View ArticleGoogle Scholar
  29. Miller KM, Li S, Kaukinen KH, Ginther N, Hammill E, Curtis JMR, et al. Genomic signatures predict migration and spawning failure in wild Canadian salmon. Science. 2011;331:214–7.View ArticlePubMedGoogle Scholar
  30. Fausch KD, Torgersen CE, Baxter CV, Li HW. Landscapes to riverscapes: bridging the gap between research and conservation of stream fishes. Bioscience. 2002;52:483–98.View ArticleGoogle Scholar
  31. Ficke AD, Myrick CA, Hansen LJ. Potential impacts of global climate change of freshwater fisheries. Rev Fish Biol Fisher. 2007;17:581–613.View ArticleGoogle Scholar
  32. Gustafson RG, Waples RS, Myers JM, Weitkamp LA, Bryant GJ, Johnson OW, et al. Pacific salmon extinctions: quantifying lost and remaining diversity. Conserv Biol. 2007;21:1009–20.View ArticlePubMedGoogle Scholar
  33. Perry AL, Low PJ, Ellis JR, Reynolds JD. Climate change and distribution shifts in marine fishes. Science. 2005;308:1912–5.View ArticlePubMedGoogle Scholar
  34. Wenger SJ, Isaak DJ, Luce C, Neville HM, Fausch KD, Dunham JB, et al. Flow regime, temperture, and biotic interactions drive differential declines of trout species under climate change. Proc Nat Acad Sci. 2011;108:14175–80.View ArticlePubMed CentralPubMedGoogle Scholar
  35. Isaak DJ, Rieman BE. Stream isotherm shifts from climate change and implications for distributions of ectothermic organisms. Glob Change Biol. 2013;19:742–51.View ArticleGoogle Scholar
  36. Allendorf FW, Hohenlohe PA, Luikart G. Genomics and the future of conservation genetics. Nat Rev. 2010;11:697–709.View ArticleGoogle Scholar
  37. Cooke SJ, Sack L, Franklin CE, Farrell AP, Beardall J, Wikelski M, et al. What is conservation physiology? Perspectives on an increasingly integrated and essential science. Conserv Physiol. 2013;1:1–23.View ArticleGoogle Scholar
  38. Narum SR, Buerkle CA, Davey JW, Miller MR, Hohenlohe PA. Genotyping-by-sequencing in ecological and conservation genomics. Mol Ecol. 2013;22:2841–7.View ArticlePubMed CentralPubMedGoogle Scholar
  39. Hoffmann AA, Sgro CM. Climate change and evolutionary adaptation. Nature. 2011;470:479–85.View ArticlePubMedGoogle Scholar
  40. Portner HO. Oxygen- and capacity-limitation of thermal tolerance: a matrix for integrating climate-related stressor effects in marine ecosystems. J Exp Biol. 2010;213:881–93.View ArticlePubMedGoogle Scholar
  41. Pandey RV, Franssen SU, Futschik A, Schlotterer C. Allelic imbalance metre (Allim), a new tool for measuring allele-specific gene expression with RNA-seq data. Mol Ecol Res. 2013;13:740–5.View ArticleGoogle Scholar
  42. Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009;10:R25.View ArticlePubMed CentralPubMedGoogle Scholar
  43. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26:139–40.View ArticlePubMed CentralPubMedGoogle Scholar
  44. Dillies M-A, Rau A, Aubert J, Hennequet-Antier C, Jeanmougin M, Servant N, et al. A comprehensive evaluation of normalization methods for Illumina high-throughput RNA sequencing data analysis. Brief Bioinform. 2012;10:1093.Google Scholar
  45. McCarthy DJ, Chen Y, Smyth GK. Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation. Nucleic Acids Res. 2012;40:4288–97.View ArticlePubMed CentralPubMedGoogle Scholar
  46. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc B. 1995;57:289–300.Google Scholar
  47. Conesa A, Götz S, Garcia-Gomez JM, Terol J, Talon M, Robles M. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics. 2005;21:3674–6.View ArticlePubMedGoogle Scholar

Copyright

© Narum and Campbell; licensee BioMed Central. 2015

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Advertisement