Transcriptomic responses to environmental change in fishes: Insights from RNA sequencing

The need to better understand how plasticity and evolution affect organismal responses to environmental variability is paramount in the face of global climate change. The potential for using RNA sequencing (RNA-seq) to study complex responses by non-model organisms to the environment is evident in a rapidly growing body of literature. This is particularly true of fishes for which research has been motivated by their ecological importance, socioeconomic value, and increased use as model species for medical and genetic research. Here, we review studies that have used RNA-seq to study transcriptomic responses to continuous abiotic variables to which fishes have likely evolved a response and that are predicted to be affected by climate change (e.g., salinity, temperature, dissolved oxygen concentration, and pH). Field and laboratory experiments demonstrate the potential for individuals to respond plastically to shortand long-term environmental stress and reveal molecular mechanisms underlying developmental and transgenerational plasticity, as well as adaptation to different environmental regimes. We discuss experimental, analytical, and conceptual issues that have arisen from this work and suggest avenues for future study.


Introduction
There are two primary mechanisms by which animal populations might alter phenotypes in response to environmental change.Plasticity (the ability of a genotype to produce different phenotypes, depending on environmental conditions; Bradshaw 1965) shifts the trait phenotype along a "norm of reaction" (sensu Woltereck 1909) defined by the genotype.Evolutionary (i.e., genetic) change occurs when selection acts on standing genetic variation to alter allele frequencies.These mechanisms are not mutually exclusive; for example, plasticity can facilitate adaptive evolutionary change to new environments (Ghalambor et al. 2007(Ghalambor et al. , 2015)), and the magnitude and direction of plastic responses can themselves evolve in response to selection (Bradshaw 1965;Schlichting 1986;Lande 2009;Chevin et al. 2010).The need to better understand how plasticity and evolution will affect the responses of wildlife to environmental variability is intensifying as the magnitude and inevitability of global climate change becomes increasingly clear (Merilä and Hendry 2014).
As with other challenges related to climate change, technological advances will likely play a key role in addressing its effect on global biodiversity.Transcriptomics (the study of all the genes expressed at a particular moment; Wang et al. 2009) has been used to unravel the relationships between environment, genotype, and phenotype in natural populations for over a decade.Hybridization-based microarrays have traditionally been the dominant method for characterizing genome-wide expression levels in ecological studies (Alvarez et al. 2015).However, RNA sequencing (RNA-seq) is being increasingly used as next-generation sequencing technologies become more accessible (refer to Box 1 for an outline of some of the key features and limitations of these methodologies in the context of this review).Although still in its relative infancy, a remarkably rapidly growing body of literature has capitalized on the advantages of RNA-seq to study how fishes (the most speciose of vertebrates) are impacted by environmental change.This work is motivated by the ecological and socioeconomic value of wild fish populations globally, the aquaculture industry, and their popularity as model systems for medical and genetic research.Yet, a unique feature of RNA-seq as it applies to ecological questions about global change is the sheer diversity of species whose investigation is unlocked through de novo transcriptome assembly.An ISI Web of Science search using the keywords "fish" and "RNA-seq*" from 2008 (when the first studies using the term "RNA-seq" were published) returned 605 articles, the majority (61%) of which were published in the previous two years.Of the 420 studies remaining after excluding those that focused on organisms other than fishes, many examined the effects of diseases (56), parasites (7), pollutants (45), and diet (19) (see Qian et al. 2014).Although these topics are relevant to the discussion of responses to environmental change because of interactions between biotic and abiotic factors, herein, we focus solely on direct responses by fishes to well-studied abiotic environmental variables; i.e., those to which animals have evolved a response, but to which changes in climate are expected to alter the conditions that they naturally experience.Although this might seem unduly restrictive, such studies proffer unique opportunities to advance our understanding of plastic and evolutionary responses to environmental change.In addition, they are often accompanied by specific challenges that we feel can be ameliorated by an early critical review of this rapidly expanding field.
After including an additional six studies not detected by the keyword search, the resultant 52 studies investigated the responses of a remarkably diverse 38 fish species to one or more of the following factors: temperature (24), salinity (20), dissolved oxygen concentration (11), and pH (7) (Table 1), all of which are expected to be affected in the coming decades as a result of global climate warming.Water temperatures are increasing across most of the globe, whereas pH and dissolved oxygen concentration are declining (IPCC 2013).Salinity is increasing in some freshwater and oceanic environments where evaporation outpaces precipitation (IPCC 2013;Settele et al. 2014) and decreasing in some tropical and high latitudes due to increased precipitation and sea ice melt (Durack et al. 2012).We review these studies in the context of how RNA-seq can inform our understanding of plastic and evolutionary responses to environmental change in fishes.We note that it is premature to draw firm conclusions based on the limited body of work thus far, especially considering variation among the tissues and time points sampled as these can drastically impact expression.However, in the interest of inspiring future investigations of this nature, we draw some tentative inferences based on comparative results across studies.We also highlight some of the experimental, analytical, and conceptual issues that have arisen in the early days of using RNA-seq to study wild populations.
Box 1. Key features and limitations of microarrays, expressed sequence tags (ESTs), and RNA-sequencing.
Traditional hybridization-based (microarrays) and sequence-based (EST) transcriptomic methods suffer from some technical limitations (reviewed by Wang et al. 2009).For microarrays, which involve hybridizing fluorescently labeled complementary DNA (cDNA) with probes affixed to a solid surface (Schena et al. 1995), these include dependence on a priori knowledge of genome sequence and a narrow quantification range that is constrained by (1) high background noise generated by nonspecific hybridization (Okoniewski and Miller 2006) and (2) fluorescent signal saturation for highly abundant transcripts.Statistical treatment of microarray data has advanced to address many of its initial limitations; therefore, microarrays might still be an appropriate and cost-effective choice for quantifying differential expression in known transcripts if a microarray for the species of interest (or one closely related) is available.Although high-throughput EST methods have largely overcome the limitations of microarrays through direct cDNA sequencing, they use outdated, time-, cost-, and labour-intensive Sanger sequencing technology.Furthermore, because the tag sequences are short, incomplete transcripts, many are unable to be annotated even if a complete reference genome is available (Costa et al. 2010), and it is not possible to distinguish between alternative splice isoforms or different alleles (Wang et al. 2009).
Advances in high-throughput sequencing technology have revolutionized transcriptomics with the advent of deep RNA-seq.A rapid and cost-effective method, RNA-seq can determine suites of genes expressed in a particular environment, including their sequences with single base-pair resolution, and their relative abundances with far greater precision than previous techniques (Wang et al. 2009).Without the need for a reference genome, RNA-seq provides a relatively accessible mode for studying the complex responses of nonmodel organisms to the environment.The ability of RNA-seq to distinguish allelic and splice variants adds another layer of valuable information for this purpose.Because it allows determination of both quantitative (i.e., differential expression levels) and qualitative (i.e., sequence) variations in gene expression, RNA-seq has the potential to enable researchers to begin to disentangle the relative contributions of transcript abundance, allelic variants, and alternative splicing to phenotypic change.Transcript sequences can be used to identify single-nucleotide polymorphisms (SNPs) within coding regions either de novo (given sufficient coverage; van Belleghem et al. 2012;Gayral et al. 2013;Lopez-Maestre et al. 2016) or with the aid of a reference genome (Piskol et al. 2013).Although the ability of de novo SNP discovery to exhaustively detect all SNPs is uncertain, those identified in transcripts are arguably likely to have a direct functional impact (Lopez-Maestre et al. 2016).Variation in these SNPs can then be inexpensively and efficiently characterized in large-scale ecological and evolutionary studies (e.g., Romiguier et al. 2014).Furthermore, mRNA sequences can inform us about putative downstream protein structure and function with no prior knowledge about particular genes.RNA-seq is not without its own technical biases (e.g., fragmentation (Sendler et al. 2011) and PCR (Aird et al. 2011) biases during sample preparation), computational limitations (e.g., difficulties of assembling and aligning short reads (Engström et al. 2013)), and logistical constraints (e.g., high costs prohibiting adequate levels of replication (Table 1)).We address some of these concerns as they pertain to this review in the discussion, in which we also recommend workflows that combine multiple methods when appropriate.We anticipate that current drawbacks of RNA-seq will be ameliorated by technological and computational advances in the near future.Nevertheless, RNA-seq can be leveraged through considered experimental design to be a promising tool to address the question of how animal populations will respond to global climate change.Note: The number of biological replicates is given with the number of individuals pooled within each replicate denoted in parentheses when applicable.Key processes involved in responses include those of focal interest to the authors, those to which the greatest number of dysregulated genes were annotated, and those that showed the highest enrichment, to a maximum of five processes.Arrows indicate the relationship between environmental and response variables, where applicable.Identical superscripts between the "Experimental comparison" and "Environmental variable/challenge" columns denote which comparisons were made with which variables when multiple options are present within a study.Asterisks denote that the challenge was conducted in combination with a bacterial infection.

Short-term, acute responses
Many experiments have investigated physiological responses to short-term (herein defined as 1-120 h), often extreme conditions, usually with the aim of discovering "tolerance" genes for breeding or aquaculture (e.g., Liu et al. 2013;Xia et al. 2013;Ao et al. 2015;Sun et al. 2015).These challenge studies reveal genes involved in short-term stress responses, many of which quickly return to baseline levels with no consequences to fitness (van Straalen and Feder 2012).These genes are arguably less likely to be targets for selection during the expected gradual, directional environmental shifts associated with climate change (Kassahn et al. 2007;Logan and Somero 2010).
Nonetheless, these studies can be informative about the way environmental signals are integrated and how response pathways evolve, and might help us to understand the resiliency of populations in the face of extreme weather events that are expected to increase in frequency and magnitude (Rahmstorf and Coumou 2011).An extensive and rapid response of the transcriptome to high salinity was found in euryhaline medaka (Oryzias latipes;  (Clark et al. 2008).These studies point towards the existence of yet another heat coping mechanism that warrants further study.As with osmoregulatory stress, the effect of heat stress on metabolic processes varied among species (Table 1).Notably, metabolism was strongly upregulated in the heat-stressed catfish (Liu et al. 2013) and delta smelt (Hypomesus transpacificus; Jeffries et al. 2016), species known for their high thermal tolerance, whereas a lack of metabolic response was observed in P. borchgrevinki (Bilyk and Cheng 2014).Such evidence of contrasting responses to heat stress in cold-and warm-adapted fishes sheds light on how adaptive divergence can alter the contents of the genomic tool kits with which species can respond to contemporary thermal stress.In the model zebrafish (Danio rerio), acute cold stress has repeatedly been associated with large shifts in transcriptional regulation (Long et al. 2013(Long et al. , 2015;;Hu et al. 2015a;Hung et al. 2016).This stress response substantially overlaps with that induced by hypoxia, both involving the upregulation of many genes involved in oxygen transport (Long et al. 2013;Long et al. 2015).Hypoxia tolerance has been associated with variation in the expression of genes involved in a variety of processes, including the regulation of epithelial permeability (Sun et al. 2015) and repression of cellular apoptosis (Yuan et al. 2016) in channel catfish (Ictalurus punctatus), avoiding cerebral inflammation in the large yellow croaker (Larimichthys crocea; Ao et al. 2015), and lipid utilization in a hybrid striped bass (Morone sp.; Beck et al. 2016).
This limited collation of studies suggests a global coordination of stress response in teleost fishes combined with the regulation of stress-specific genes dependent on species-specific adaptations.

Long-term, chronic responses
After the initial stress response, how do fishes adjust physiologically (i.e., acclimatize) during prolonged exposure to new environmental conditions?Given that climate change will involve sustained alteration of the environment, the genes and pathways identified in long-term experiments (herein defined as 1-4 weeks) are more likely to be involved in some form of a plastic response that has fitness consequences (either adaptive or maladaptive) (Smith et al. 2013).
The general stress response is less apparent following prolonged exposure to increased salinity, consistent with the hypothesis that acclimation is common (Table 1).A vast array of genes and pathways has been proposed to enable prolonged salinity tolerance, including those involved in ion transport (Lam et al. 2014;Wang et al. 2014;Gu et al. 2015;Nguyen et al. 2016), blood pressure regulation and fat metabolism (Xu et al. 2015), and both innate and adaptive immunity (Norman et al. 2014a(Norman et al. , 2014b)).
The mechanisms by which long-term hypoxia leads to reproductive impairment in the marine medaka (Oryzias melastigma) have been determined through sex-specific brain transcriptome sequencing (Lai et al. 2016b) coupled with gonadal microRNA profiling (Lai et al. 2016a;Tse et al. 2016).Hypoxia-responsive microRNAs (small non-coding RNAs which can post-transcriptionally modulate gene expression; Carrington and Ambros 2003) were associated with the upregulation of steroidogenic enzymes and hormone receptors in the ovary (Lai et al. 2016a) and diverse cellular processes including epigenetic modifications in the testes (Tse et al. 2016).
HSPs and immune-related genes associated with the short-term heat stress response are likewise upregulated during prolonged heat exposure in crimson-spotted rainbowfish (Melanotaenia duboulayi; Smith et al. 2013), redband trout (Oncorhynchus mykiss gairdneri; Narum and Campbell 2015), and half-smooth tongue sole (Cynoglossus semilaevis; Guo et al. 2016) whereas metabolic processes continue to be one of the most enriched categories of dysregulated genes during long-term heat stress in these studies.However, immunity-related genes comprise a much smaller proportion of differentially expressed transcripts in rainbowfish when compared with short-term challenge studies in other species.Further, both the number of differentially expressed transcripts and the expression levels of stress response genes decreased over the 4 week duration of the study on redband trout (Narum and Campbell 2015), suggesting acclimation to heat stress.Limited evidence from microarrays cautions that even if the stress response decreases following acclimation in heat-tolerant species, maintenance costs for homeostasis might be higher at warmer temperatures (Logan and Buckley 2015): energetically costly protein biosynthesis and active ion transport were upregulated in the longjaw mudsucker (Gillicthys mirabilis) after three weeks of heat exposure, whereas HSPs were largely absent (Logan and Somero 2010).Less energy for foraging, growth, and reproduction would be available to species with such a response.If similar biological processes characterize the majority of dysregulated pathways in both acute and long-term responses to environmental change, then challenge experiments might reliably be used to uncover the general physiological processes underlying long-term plastic responses.However, there is evidently potential for acclimation to decrease the magnitude of the plastic response, and it is not clear whether the same genes are involved at different times during the response.

Developmental plasticity
Environmental conditions experienced earlier in life can both alter subsequent phenotypes and impact future plastic responses to the environment through epigenetic mechanisms (i.e., those which "cause chromosome-bound, heritable changes to gene expression that are not dependent on changes to DNA sequence"; Deans and Maggert 2015, p. 892).This developmental plasticity or acclimation can enhance persistence in new and variable environments and result in novel phenotypes that can facilitate adaptation (West-Eberhard 2003).
For example, through RNA-seq, we are beginning to better understand the mechanisms underlying the ability of thermal acclimatization to shift the breadth and optima of thermal performance later in life (sensu Fry and Hart 1948).Embryonic exposure to thermal extremes appears to enhance the response of adult zebrafish to cold temperatures, surprisingly resulting in greater swimming performance regardless of the direction of the extreme (Scott and Johnston 2012).The improved acclimation capacity of the warm-incubated fish was explained by differential expression of genes involved in energy metabolism, blood vessel development, and muscle contraction and remodelling, which corresponded with differences in muscle area and composition (Fig. 1).It would be interesting to know whether the cold-incubated fish (which were not sequenced) achieved acclimation via the same transcriptional modifications.In a separate study, RNA-seq revealed the molecular basis of "rapid cold hardening" (Kelty and Lee 2001), whereby brief exposure to mild cold improved larval survival in the face of severe cold stress (Long et al. 2013).Promoter switching and alternative splicing emerged as major mechanisms enabling cold tolerance in fishes, consistent with previous studies on a wide range of stressors in other taxa, although the functional significance of different isoforms remains to be investigated.

Transgenerational plasticity
Non-genetic parental influences on offspring phenotype can facilitate acclimation across generations (Mousseau and Fox 1998).Evidence of such transgenerational effects in fishes suggests that they might play a major role in enabling fish populations to cope with environmental change (Donelson et al. 2012;Hurst et al. 2012;Miller et al. 2012;Salinas and Munch 2012), particularly in species that have less capacity for acclimation as adults because they have evolved in a relatively stable environment (e.g., coral reef fishes; Munday et al. 2012).Veilleux et al. (2015) explored the molecular basis of this phenomenon, using RNA-seq in a common reef fish (Acanthochromis polyacanthus), by evaluating gene expression and metabolic performance in response to increased temperature both within and across generations.Differential expression was greater in transgenerationally exposed fish, which had improved aerobic scope, compared with developmentally exposed fish, for which performance was reduced relative to the controls (Fig. 2).The biological processes associated with the developmental response to temperature were also part of the transgenerational response (e.g., metabolism, immunity, and stress response), as were genes previously found to respond to short-term thermal challenge (e.g., apolipoproteins; Kassahn et al. 2007;Podrabsky and Somero 2004), suggesting a link between short-term, developmental, and transgenerational thermal stress responses in fishes.Interestingly, HSPs were largely absent from both developmental and transgenerational treatments, suggesting that they might not be good predictors of thermal acclimation capacity.The epigenetic mechanisms involved in regulating the developmental and transgenerational thermal responses above are unknown.In another study, microRNA sequencing has revealed a specific epigenetic effect of hypoxia that causes transgenerational reproductive impairments in male marine medaka (Tse et al. 2016).Sequencing of microRNAs is, therefore, another promising avenue for understanding how gene expression is fine tuned by epigenetic mechanisms in response to environmental factors throughout development and across generations.The epigenetic mechanisms responsible for regulating developmental and transgenerational plasticity are of substantial interest given their considerable potential to improve our understanding of the capacity of fishes to cope with rapid environmental change.

Responses to multiple stressors
Climate change is altering many environmental variables simultaneously (IPCC 2013).Considering that multiple stressors can have complex interactive effects (Schulte 2007), studies examining the combined effects of heat and other stressors on fishes are highly relevant to predictions of fish responses to global climate change.
Heat stress suppressed the immune system of Mozambique tilapia (Oreochromis mossambicus) infected with a bacterial pathogen, apparently through metabolic constraints imposed by limited oxygen (Wang et al. 2016).Among the 2000+ differentially expressed genes, rates of synonymous and nonsynonymous substitutions based on SNPs identified from the O. mossambicus transcriptome and the closely related, but less disease resistant, Nile tilapia (Oreochromis niloticus) revealed signs of positive selection in O. mossambicus for 43 genes involved in the immune response and oxidative respiration.These findings suggest that O. mossambicus has evolved superior disease resistance relative to O. niloticus, yet its immune system is impaired by heat stress.A better understanding of how temperature mediates infection in fishes, many of which have unusual or poorly understood immunological strategies (Buonocore and Gerdol 2016), is urgent as climate change increases the incidence of disease outbreaks globally (Brander 2007).
Along with rising temperatures, acidification driven by increases in dissolved carbon dioxide is a major threat to fishes (Pörtner et al. 2004).A long-term dual-stressor time-series experiment on the Antarctic notothenioid P. borchgrevinki suggests that, when occurring in tandem, these shifts can produce distinct responses when compared with heat stress alone (Huth and Place 2016a).Huth andPlace (2016a, 2016b) demonstrated an inflammatory response to increased temperatures and pCO 2 that lasted at least 7 d, along with an increase in rates of cell death followed by gradual acclimation to near basal expression levels by 56 d, in two Antarctic notothenioids.However, the degree of response was reduced overall in P. borchgrevinki compared with T. bernacchii, suggesting that sensitivity to environmental perturbation varies among these closely related cold specialists.In contrast, the long-term response of Amazonian tambaqui to these dual stressors was dominated by molecular chaperones and metabolic and developmental processes (Prado-Lima and Val 2016).
The interaction between multiple abiotic stressors has been explored using RNA-seq from a developmental perspective with respect to the extent to which variation in prior exposure along one environmental axis influences the response along a different axis.Although developmental cold exposure has been reported to protect larval zebrafish against future cold stress, it was associated with decreased tolerance to lethal hypoxia, whereas prior exposure to mild hypoxia improved both hypoxia and cold tolerance (Long et al. 2015).Genes involved in oxygen transport were mainly associated with this process, revealing molecular mechanisms underlying the hypothesis that oxygen limitation is the primary determinant of thermal tolerance in fishes (Pörtner 2002).Somewhat less intuitively, acclimation to different salinities activated different strategies to cope with cold tolerance in milkfishes (Chanos chanos), whereby seawater-acclimated milkfish were more cold tolerant than those acclimated to freshwater (Hu et al. 2015b).The seawater-acclimated fish upregulated a suite of genes related to increasing the energy budget, whereas freshwater-acclimated fish reduced energy loss by downregulating genes involved in basal metabolism.
These studies highlight the fact that previous exposure and interactions between multiple stressors can have substantial and perhaps surprising consequences on fitness and are thus critical to understanding fish responses to our changing climate.

Evolutionary responses to environmental change
Although individual and transgenerational plasticity can help organisms cope with environmental change in the short term (e.g., fewer than five generations), responding to ongoing climatic shifts will involve evolution because traits are not necessarily plastic and (or) reaction norms will no longer be adaptive in the new environment (Visser 2008).Although empirical evidence is scarce relative to that of plastic responses (Merilä and Hendry 2014), rapid evolution in response to environmental change has been documented in a variety of taxa (e.g., Bradshaw and Holzapfel 2001;Umina et al. 2005;Derry and Arnott 2007;Charmantier et al. 2008), including fishes (reviewed by Crozier and Hutchings 2014).This section summarizes what has been learned about adaptive responses to environmental change in fishes from RNA-seq experiments.

Identifying candidate genes for adaptation
A primary aim of transcriptomics is to identify candidate genes for adaptation; i.e., those genes with large impacts on fitness under different environmental conditions (Feder and Mitchell-Olds 2003).This "discovery-driven" approach proved powerful early on in the study of genomic reaction norms (Aubin-Horth and Renn 2009).RNA-seq offers an advantage for candidate gene discovery because of its unbiased nature and lack of necessity for prior information (Wang et al. 2009).The studies described previously identified numerous candidate genes that are potential targets for selection in response to changes in temperature, salinity, dissolved oxygen concentration, and pH.Armed with such information, researchers can develop functional markers to monitor for a contemporary response to climate change (Hoffmann and Willi 2008) or screen broadly across a species range to predict the potential for adaptation (Hoffmann and Sgrò 2011).The unique opportunities proffered by RNA-seq have yet to be fully taken advantage of with regard to the evolutionary effects of environmental change on fishes, but this avenue of research holds great potential.

Intraspecific variation in transcriptomes
Transcriptomic variation at the population level can reveal how gene expression evolves in response to local environmental regimes.Zhang et al. (2015) compared the transcriptomes of two ecotypes of scaleless carp (Gymnocypris przewalskii) from saline and freshwater lakes.Of the many thousands of genes that were differentially expressed, they used sequence information to narrow in on just 242 protein-coding genes that showed signs of strong positive selection.The authors concluded that relatively few genes, chiefly those involved in ion regulation and the immune response, play critical roles in the shift from saline to freshwater habitats in fishes.An acclimation experiment comparing native freshwater to anadromous saltwater threespine stickleback (Gasterosteus aculeatus) also revealed many genes potentially underlying salinity adaptation (Wang et al. 2014).Finally, a rare study examining the effects of elevated pH on a fish transcriptome described how changes in gene expression played a key role in the relatively recent shift of Amur ide (Leuciscus waleckii) from freshwater to extreme alkalinity in a soda lake (Xu et al. 2013).Further research is needed to determine whether the intraspecific variation in expression described for the wild-caught fish used in these studies has a genetic (as opposed to epigenetic) basis, thereby representing adaptive evolution.This could be achieved through traditional labour-intensive common-garden experiments or some genetic inferences could perhaps be gleaned more readily by examining allele-specific expression patterns.

Population-level variation in transcriptional plasticity
In addition to shifting the mean phenotype (i.e., reaction norm elevation), evolutionary responses to environmental change can alter the shape of the plastic response (i.e., reaction norm slope) (Bradshaw 1965;Lande 2009;Chevin et al. 2010).In fishes, there is evidence from common-garden experiments of population-and family-level variation in plastic responses to each of the environmental variables discussed thus far (reviewed by Hutchings 2011 andOomen andHutchings 2015).These experiments are extremely useful for detecting adaptation when combined with measures of fitness; however, they do not tell us about the genetic mechanisms underlying plastic responses.
RNA-seq allows us to bridge the gap between genotype and phenotype by linking genetic variation directly to differences in gene expression and then to phenotypic responses observed in the lab.Such a complete chain has yet to be made within the scope of the present review.Narum and Campbell (2015) came the closest when they found population-specific patterns of plasticity in response to heat stress among desert and montane redband trout.The desert population exhibited greater differential expression (in both the number of genes and the magnitude of the fold change) compared with either the montane population or their F1 cross (Fig. 3), although all populations showed evidence of acclimation during the 28-d experiment.HSPs were not upregulated as much in the desert trout when compared with the montane trout, whereas many genes involved in metabolic and cellular processes were highly upregulated, suggesting that the desert trout have evolved complex and specialized molecular mechanisms to cope with heat stress.The F1 cross generally exhibited intermediate expression patterns between the two populations, consistent with additive genetic variation, although a greater number of shared differentially expressed genes with the maternal montane trout suggests a possible maternal or dominant effect at some genes.
Although this novel experiment provides considerable insight into the molecular basis of thermal adaptation, corresponding physiological and other phenotypic measurements would be extremely valuable in understanding how changes in individual gene expression are related to fitness.Further, the inclusion of multiple temperatures in such an experiment would allow for direct quantification of thermal plasticity within populations.This approach was used by Morris et al. (2014) to compare levels of transcriptomic thermal plasticity, as determined by microarrays, between ancestral marine and derived freshwater threespine stickleback populations.More plastic genes were detected in the derived freshwater populations, supporting the hypothesis that if greater environmental variability is encountered following colonization of new habitats, it will drive the evolution of greater plasticity relative to the ancestral population.

Family-level variation in transcriptional plasticity
In addition to population-level variation in transcriptional plasticity, studies in fishes have shown differences in patterns of gene expression at the family level.Norman et al. (2014b) used a quantitative trait locus (QTL) approach to reveal the genetic basis underlying the correlation between salinity tolerance and differential expression of immune genes among families of Arctic char (Salvelinus alpinus).Interestingly, the majority of QTL associated with ion transport were located near differentially expressed genes, suggesting that cis-regulatory elements (non-coding DNA that regulates transcription of nearby genes) are involved in controlling their expression.Concomitantly, the majority of differentially expressed genes were not associated with QTL, suggesting that they might be controlled by trans-regulatory elements (those on distant genes).That genetically based differences in gene expression were found within a lab-bred strain suggests that there might be substantial variation in transcription within populations.Knowledge of the extent to which such transcriptional variation is heritable is thus far rather limited, as attempts to quantify heritability of gene expression are rare.Critically, the ability of a population or species to respond evolutionarily to environmental change depends on the amount of heritable genetic variation it possesses for adaptive traits (Hoffmann and Sgrò 2011).In a unique example, McCairns et al. (2016) used crimson-spotted rainbowfish of known pedigree to demonstrate moderate levels of heritability for transcription and transcriptional plasticity in thermal responses of candidate genes that were previously identified using RNA-seq (Smith et al. 2013).Abundant family-level variation in plasticity for most genes suggests that substantial heritable variation in plasticity might exist within fish populations for selection to act upon to produce an adaptive evolutionary response to climate change.

Challenges and directions for future research
The application of RNA-seq to study non-model systems is still in its infancy, and care must be taken not to overlook its limitations in the rush to adopt the technology (Costa et al. 2013;Todd et al. 2016).In this section, we describe some of the experimental, analytical, and conceptual challenges raised specifically by the body of work described herein and highlight some avenues of future exploratory and confirmatory analyses.For further detail about technical and analytical issues surrounding RNA-seq experiments in general, we refer the reader to several discussions on the topic (Fang and Cui 2011;Ozsolak and Milos 2011;Sendler et al. 2011;Costa et al. 2013;Conesa et al. 2016).

Experimental protocols and sampling design
Gene expression is extremely sensitive in that it can potentially change very quickly and in response to slight environmental variability.This can make sampling protocols challenging, as any manipulations that fish experience between the environment of interest (whether in the field or laboratory) and the preservation of tissue for RNA extraction can affect gene expression.Depending on the research question, it might be prudent to collect a control sample for handling or transfer stress (e.g., Wong et al. 2014) or avoid feeding prior to sampling (e.g., Liu et al. 2013).
Although relatively inexpensive from an input:output (cost:amount of data generated) perspective, the substantive financial, laboratory, and computational resources required for RNA-seq appear to have limited experimental design in terms of sample sizes and numbers of replicates and treatments.Indeed, of the studies included here, nearly half did not include true biological replicates (at best, samples were pooled within groups to create a single library for RNA-seq), consistent with the dearth of replication observed in ecological and evolutionary studies in general (Todd et al. 2016).St. Laurent et al. (2013) concluded that as few as three biological replicates can be sufficient to detect small differences in expression in a model system, whereas Robles et al. (2012) found that at least six are necessary when dealing with low count data (e.g., when using a multiplex sequencing strategy).However, the increased biological variation inherent in studies on non-model organisms essentially means that the more biological replication, the better, as it is the only way to cope with the unique statistical challenges of controlling for false positives during such highly parallelized significance testing while retaining the ability to detect lowly expressed transcripts.Although few replicates can still produce significant results, they will be biased towards genes with large log-fold changes, and these are not necessarily the most relevant for fitness (Smith et al. 2013;Evans 2015).
Few of the experimental RNA-seq studies published to date included biological replicates within tank replicates, which is necessary for estimating variation due to random tank effects.This might be partly explained by the difficulty of such computations with currently available differential expression analysis software.However, a few experiments dispersed replicates among tanks to eliminate tank variation as a confounding factor (e.g., Narum and Campbell 2015;Prado-Lima and Val 2016).Studies measuring multiple time points are also uncommon because they require multiplying sample sizes, yet they are key to distinguishing between short-term stress and long-term acclimation responses (e.g., Narum and Campbell 2015;Xu et al. 2015;Huth andPlace 2016a, 2016b).The number of individual and tank replicates required is highly dependent on the amount of biological variation within and between groups, which can be difficult to anticipate as it will likely be extremely case specific.However, the development of tools to help predict the power of various experimental designs based on data from pilot experiments (e.g., Scotty (Busby et al. 2013); reviewed by Todd et al. 2016), and future decreases in the cost of RNA-seq will likely facilitate more robust and complex experiments.Alternatively, the utility of RNA-seq to explore uncharacterized genomes can be leveraged to design custom microarrays (Alvarez et al. 2015) for expansion of experimental designs or to reduce the bioinformatic challenges associated with follow-up experiments.Kusakabe et al. (2017) took advantage of both RNA-seq and microarrays, combined with QTL mapping, qPCR, whole genome sequencing, and a genome scan, to confidently identify candidate genes for salinity adaptation in threespine stickleback.Notably, the direction and magnitude of expression differences for the genes identified by transcriptome analysis were similar between the RNA-seq and microarray results.

Bioinformatic analysis
The technological limitations of RNA-seq are not of as great concern when compared with the analytical limitations, as our ability to generate data vastly surpasses our ability to interpret it.Key areas for improvement in this regard are bioinformatic and genomic resources.Differential expression analysis programs capable of handling mixed-effects models are needed to control for variation due to random (e.g., tank) effects.Biological interpretation of large differentially expressed gene lists remains a major obstacle, although gene set enrichment analyses, which groups genes based on common biological function, location, or regulation (Subramanian et al. 2005), have become a popular aid in their distillation to major functional categories (e.g., Liu et al. 2013;Long et al. 2013;Gu et al. 2015).However, there is often a focus on only the largest or most expected functional categories, potentially ignoring important biological processes involving fewer genes and providing opportunity for reporting bias and subjective interpretation of results.
When the goal is to identify a reduced set of the most promising candidate genes, complementary methods using high-resolution sequence data such as testing for positive selection through comparing rates of synonymous and nonsynonymous substitutions can be helpful (e.g., Kavembe et al. 2015;Zhang et al. 2015;Wang et al. 2016).Although underutilized in ecological genomics, gene network construction offers a promising alternative for functional analysis based on the notion that genes with a higher degree of connectedness have a greater impact on fitness through involvement in multiple cellular pathways (i.e., hub genes; Costanzo et al. 2010).Although incorporation of differential expression data improves these networks by taking co-expression into account (Fu et al. 2014), it is not necessary because they are built using databases of known gene interactions (Evans 2015).If specific (i.e., candidate) genes are the unit of interest, validation of gene expression using additional samples (either more biological replicates from the same experiment or, ideally, samples from a different population or incidence of environmental exposure) would be an asset to reduce the potential for false positives.
A lack of reference data for functional annotation is one of the biggest challenges facing RNA-seq analysis in an ecological context, a consequence of the fact that annotations are based on homologies with genes from a few model species and environmental conditions (Pavey et al. 2012).For example, over half of the transcripts differentially expressed in the heat stress response of crimson-spotted rainbowfish could not be annotated (Smith et al. 2013).Pavey et al. (2012) proposed alleviating this problem by creating a database of ecological annotations (sensu Landry and Aubin-Horth 2007), allowing for cross-referencing genes and associated environmental variables across studies and species.Such a database could be populated with results from functional assays (i.e., experiments which demonstrate that a gene has a measureable effect on a particular phenotype) such as targeted gene knockdown or overexpression, as these become available.For example, Hu et al. (2015a) confirmed the coldresponsive functions of novel cis-regulatory elements discovered using RNA-seq by exposing mutant and wild-type transgenic zebrafish embryos to different temperatures.Ecological annotations would complement the current practice of using the Gene Ontology (geneontology.org) to annotate genes according to cellular component, molecular function, and biological process.However, genes associated experimentally with particular environmental variables need not necessarily be annotated to serve as functional markers for ecological studies of environmental adaptation.

Conceptual challenges
Most criticisms of RNA-seq can justifiably pertain to transcriptomics in general.For example, inferring physiological function from gene expression data is inherently problematic due to transcriptional noise (caused by nonspecific initiation of transcription by RNA polymerase II; Struhl 2007) and time lags between environmental exposure, transcriptional response, and physiological and fitness consequences (for a detailed discussion, see van Straalen and Feder 2012).With respect to its utility in discovering candidate genes of major consequence for environmental adaptation, the main criticisms (as proposed by Feder and Walser 2005) are: (1) that such genes are rare, and therefore the wide net approach of transcriptomics may not be the most effective; (2) that differences in gene expression levels do not necessarily correlate with differences in fitness; and (3) that post-transcriptional and post-translational modifications mean that mRNA abundance is not a good predictor of protein abundance, which is more relevant for fitness and can also be examined directly to address questions of environmental plasticity and adaptation (see Papakostas et al. 2014 andMäkinen et al. 2016 for examples in the European grayling (Thymallus thymallus)).Although these issues have persisted to some extent for a decade, ongoing developments from systems-level experiments, advancing bioinformatic techniques, and additional gene functional analyses will continue to improve prospects (reviewed by Evans 2015).
Specifically, increasing our ability to detect lowly expressed genes, as described in the previous section, and dedicating more effort towards interpreting these and moderately differentially expressed transcripts will advance efforts to better understand the genetic basis of plastic and evolutionary responses to environmental change (e.g., Smith et al. 2013), as plasticity need not be great to have fitness consequences.
Experiments that directly link gene expression with physiological and other phenotypic traits, particularly measures of fitness, are of vital importance for understanding how transcriptomic variation ultimately affects populations and species.One approach involves following individuals after nondestructive RNA sampling (potentially at many points during development) through to reproduction or mortality, thereby providing a more comprehensive and concrete view of how the environment interacts with development to alter gene expression and, ultimately, fitness.For example, Evans et al. (2011) combined non-lethal sampling of sockeye salmon (Oncorhynchus nerka) gill tissue and microarray expression profiling with telemetry to identify transcriptional signatures of mortality during subsequent upstream migration.In a related experiment, microarray expression profiles from non-lethal samples were compared between time-matched moribund and surviving wild-caught salmon after being held for one week at warm, ecologically relevant temperatures (Jeffries et al. 2012).Such approaches also address a caveat of most experimental RNA-based studies in which mortality occurs: that only the survivors are sampled.
Another issue inherent in "omics" approaches in general is the difficulty of testing the hypothesis that a putative candidate gene is indeed of major consequence for environmental adaptation (Evans 2015).In addition to functional validation, candidate genes derived from laboratory studies could be validated in the field such as through broad-scale screening of wild populations coupled with a priori predictions regarding the expected distributions of variants based on environmental data.A conceptually similar approach was employed to identify those RNA-seq-derived candidate genes for thermal adaptation that occur in the same molecular pathways as genes showing signs of climate-mediated selection in the wild (McCairns et al. 2016).Notwithstanding the difficulty of controlling for background genetic variation in such studies on wild populations, although feasible in some systems (e.g., ancestral versus derived populations; Barrett et al. 2008), they can independently corroborate the role of candidate genes in adaptation and contribute to a mechanistic understanding of the impact of genotype on fitness (sensu Dalziel et al. 2009).

Conclusion
The popularity of using RNA-seq to study environmental responses in natural populations has increased rapidly over the past few years and represents a promising method, alone or with other transcription quantification techniques as part of a "unifying workflow" (Alvarez et al. 2015), to understand the effects of environmental change on global biodiversity.Field and laboratory experiments on both model and non-model fishes have already provided much insight into the potential for individuals to respond plastically to short-and long-term environmental stress and for populations and species to evolve in the face of shifting environmental regimes.Numerous candidate genes for environmental adaptation have been identified for further study and myriad fish species not mentioned here have recently had their transcriptomes characterized (e.g., Li et al. 2015a;Salem et al. 2015;Salisbury et al. 2015;Yue et al. 2015;Fan et al. 2016;Kolder et al. 2016).Coupling this potential with ongoing technological and bioinformatic advances will lead to rapid developments in this field in the coming years.In particular, we expect an even greater broadening of the taxonomical and geographical representation of non-model study species and an increase in downstream analyses that use the sequence information to further test specific hypotheses derived from initial exploratory approaches.We also hope to see more studies being interpreted in an ecological or evolutionary context, as this would greatly facilitate their application to the global climate crisis.

FACETS | 2017 Fig. 1 .
Fig.1.Effect of embryonic temperature (solid line = 27 °C; dashed line = 32 °C) and long-term acclimation temperature on adult zebrafish (Danio rerio) (a) swimming performance, (b) muscle phenotype, (c) and (d) primary transcriptional responses as identified by principal component analysis, and (e) and (f) transcription of genes representative of those involved in the primary transcriptional responses, given in normalized read counts (means ± SEM).Some error bars are too small to be seen (redrawn from Scott and Johnston 2012).

FACETSFig. 3 .
Fig.3.Differential expression of (a) a desert strain, (b) an F1 cross, and (c) a montane strain of redband trout (Oncorhynchus mykiss gairdneri) exposed to heat stress versus those held at control temperatures.Significant differentiation (FDR ≤ 0.05) is indicated in red.Green and blue lines represent greater than or equal to twofold and greater than or equal to fourfold changes, respectively (original source:Narum and Campbell 2015).

Table 1 .
Research in which RNA sequencing was used to study the effect of continuous, abiotic environmental variables on fishes.

Table 1 .
(concluded ) Salinas and Munch 2012)www.facetsjournal.combyUNIVERSITETSBIBLIOTEKETi Agder on 03/15/18 later in development, or even across multiple generations (i.e., transgenerational plasticity;Salinas and Munch 2012).In this section, we discuss both short-and long-term within-generation responses, as well as non-genetic change that can occur across generations.