Effects of high temperatures on threatened estuarine fishes during periods of extreme drought

Increases in water temperatures as a result of climate change have already impacted fish populations (Perry et al., 2005; Graham and Harrod, 2009). Future climate change may substantially affect migratory fish species as direct and indirect effects of changing environmental conditions are integrated across their different life stages (Dudgeon et al., 2006; Reist et al., 2006; Crozier et al., 2008; Robinson et al., 2009). For example, increases in water temperature during development may alter the timing of seaward migration in anadromous species, causing individuals to miss optimal feeding conditions in the marine environment (Taylor, 2008). While the majority of knowledge on the effects of climate change on migratory fishes is from studying economically important anadromous species, such as Pacific salmon and Atlantic salmon (e.g. Jonsson and Jonsson, 2009; Martins et al., 2012), much less is known about the effects on anadromous and semi-anadromous species that are less important economically, yet contribute greatly to aquatic biodiversity. For anadromous and semi-anadromous fishes, estuaries provide critical rearing grounds and migratory corridors that connect the marine environment with freshwater watersheds (Levy and Northcote, 1982; Ray, 2005). Given the importance of estuaries worldwide for resident and migratory fishes, it is crucial to understand the effects of water temperature on these systems.

Estuaries can be influenced by various global and local environmental factors, in addition to anthropogenic activities and natural processes throughout the upstream watershed (Dudgeon et al., 2006; Cloern et al., 2011). Therefore, fishes reared in estuaries may be impacted by multiple environmental stressors that include habitat loss, nutrient and contaminant inputs, and colonization of invasive species (Kemp et al., 2005; Brooks et al., 2012; Cloern and Jassby, 2012). These stressors can interact with increases in water temperature, making it difficult to predict the consequences of climate change on estuarine fishes. The San Francisco Estuary/Sacramento–San Joaquin Rivers Delta system, California, USA, hereafter called the ‘Delta’, is a heavily altered estuary that has been impacted by the aforementioned stressors, in addition to substantial water diversions from the system for human use (Sommer et al., 2007). At the crossroads of human demands for water and conservation of the Delta are several high-profile native pelagic fishes (Sommer et al., 2007). These include two estuary-dependent Osmeridae fishes, the threatened anadromous longfin smelt (Spirinchus thaleichthys), the southernmost population of the species that ranges from California to Alaska, and the endangered delta smelt (Hypomesus transpacificus), a semi-anadromous species with a distribution confined to the Delta (Moyle, 2002; Sommer and Mejia, 2013). Abundances have reached record lows in California, and conservation efforts to protect these species by regulating water withdrawals from the system are particularly contentious as this affects the water supply for >25 million people and a multibillion dollar agriculture industry (Brown et al., 2013). California is in the midst of a historic drought [2012–2016 (ongoing)] that has further reduced freshwater availability (Diffenbaugh et al., 2015). California normally experiences intense wet or dry periods; however, 2014–2015 was exceptionally dry, with the drought status reaching its most severe level, leading to high surface water temperatures (Fig. 1). It is feared that extreme drought may exacerbate the declines of longfin smelt and delta smelt and potentially lead to extirpations of these species from the wild (Peter B. Moyle, University of California, Davis, personal communication).

Climate change may intensify the declines of endangered and threatened species, especially those with limited distributions, impacting aquatic biodiversity (McMahon and Hays, 2006; Schwartz et al., 2006; Carpenter et al., 2008). The body temperatures of aquatic ecotherms conform to ambient water temperatures (Beitinger et al., 2000), and species living near their thermal limits may be particularly sensitive to future climate warming (Somero, 2010). Therefore, understanding the thermal limits and physiological consequences of elevated water temperature is crucial for determining the impacts of climate change on sensitive species. However, studying the physiological responses of endangered and threatened species is challenging because specimens needed to obtain sufficient samples are often limited (McMahon and Hays, 2006). In some instances, researchers may have the opportunity to study imperiled fishes derived from wild-collected stock and reared in captivity, which allows for the examination of physiological responses to environmental stressors directly on sensitive, endangered fishes (Komoroske et al., 2015). This can be preferable to extrapolating responses of surrogate species to those of immediate conservation concern, as many stress responses can be species-specific (e.g. Jeffries et al., 2015). Understanding how high water temperatures affect threatened and endangered estuarine fishes is a crucial research focus with major ramifications for global species conservation.

Increasingly powerful approaches in the study of non-model species have come from RNA-sequencing technology (Ekblom and Galindo, 2011). Consequently, RNA sequencing has been used to characterize the effects of high water temperatures on various ecologically relevant non-model fishes (Smith et al., 2013; Narum and Campbell, 2015; Veilleux et al., 2015). Transcriptome-wide screening approaches provide valuable insight into the complex responses to water temperature, characterizing general responses to temperature in fishes (i.e. heat shock response; Narum and Campbell, 2015) and the potentially regional or species-specific responses (i.e. immune responses; Jeffries et al., 2012a). Ideally, these cellular responses can be linked with whole-organism and population-level performance indices to aid in predicting the longer-term consequences of high water temperatures on fishes. Because no prior molecular information is required for transcriptome sequencing, de novo transcriptome assembly approaches may be an extremely valuable tool for studying the effects of environmental stressors on endangered or threatened fishes that have been relatively understudied from a physiological perspective, or that are not amenable to standard experimental protocols.

As climate change scenarios will likely result in reduced precipitation along with an increase in water demand from a growing population, high water temperatures associated with drought and low water levels may represent ‘normal’ conditions in the future (Cloern et al., 2011). In the present study, we used an integrative approach to assess the effects of elevated water temperatures on the longfin smelt, a species with no prior information available on the physiological responses to water temperature, and compared these responses with those of the sympatric delta smelt. We used RNA sequencing to characterize the cellular responses to temperature in both species. We also examined differences in oxygen consumption and whole-organism thermal tolerances to link transcriptomic changes to functional responses in whole-organism performance. Such approaches may further the understanding of the functional significance of changes in transcription profiles (e.g. Nikinmaa et al., 2013). These physiological responses were coupled with temperature distributions in the wild to assess the relative sensitivity of these species to environmentally relevant high water temperatures. Based on life history differences and catch data, we hypothesized that longfin smelt would have a lower thermal tolerance than delta smelt and that responses would be detectable across multiple levels of biological organization (i.e. cellular to whole-organism responses). The temperature exposures were conducted on juveniles as developing fishes may be especially sensitive to elevated water temperatures (Johnston, 2006). To our knowledge, this is the first comparative genomics study on threatened and endangered estuarine fishes, and provides a comprehensive data set on the potential effects of climate change on two fishes at the forefront of balancing ecosystem needs and human water demand, relevant to the conservation of estuarine fishes and estuaries globally.

To describe the effects of drought on water temperatures and to characterize the thermal distributions for delta smelt and longfin smelt in the Delta, water temperatures and catch data from 34 sampling sites were obtained from the California Department of Fish and Wildlife (https://www.wildlife.ca.gov). Wild collection data are from the ‘20 mm Survey’ conducted between March and July from 1995 to 2015 (Sommer and Mejia, 2013), which collects the life-stage most relevant for this study. The delta smelt collected by the 20 mm survey were 21.2 mm (±9.58 mm s.d.; n=12,656) and the longfin smelt were 18.5 mm (±6.29 mm s.d.; n=23,382) in length. There can be significant saltwater intrusion into the Delta; however, these sites were predominantly freshwater based on salinity measurements during data collection. The water temperatures from 1995 to 2015 were collected at the same locations as the fish. We compared the March–July temperature trends, the time of year most relevant for this life stage, between 1995–2011 and 2012–2015 (drought years) to assess the effects of extreme drought on the habitat of larval and juvenile delta smelt and longfin smelt.

Wild, maturing longfin smelt were collected from the Delta by the United States Fish and Wildlife Service in December 2013 (temperature exposure studies) and 2014 (metabolic rate trials). Fish were transported to the University of California, Davis Fish Conservation and Culture Laboratory and reared until reproductively mature. These adults were used to produce the larval fish for the present study. Maturing F6 delta smelt were spawned in 2013 from a captive population maintained by the culture facility to produce the larval delta smelt used in the present study. The delta smelt breeding program incorporates a genetic management strategy to minimize inbreeding, maintain genetic representation of the wild population and maximize genetic diversity (Fisch et al., 2013). Longfin smelt and delta smelt were reared at 12°C and 16°C, respectively, until they were ∼40 days post-hatch (dph), when they were transferred to 130 liter acclimation tanks that were gradually brought to 14°C over 2 days. Age-matched longfin smelt and delta smelt were then held at 14°C for 8 days prior to the temperature exposures. A subset of fish were measured to estimate the size range of individuals used in the study (longfin smelt, n=17, standard length=15.0–19.0 mm, mass=0.006–0.014 g; delta smelt, n=43, standard length=14.8–22.0 mm, mass=0.010–0.038 g). This research was approved by the University of California, Davis Institutional Animal Care and Use Committee (protocols 17522 and 18175).

We determined the upper acute temperature tolerance of 50 dph longfin smelt (n=17) and delta smelt (n=20) using critical thermal maximum (CT_max) methodology (detailed methods in Komoroske et al., 2014). Briefly, for each CT_max trial, an individual fish was placed into a 1 liter glass container at the acclimation temperature of 14°C. The water temperature was raised at a rate of 0.3°C min⁻¹ until there was loss of equilibrium, and this temperature was used as the CT_max for that individual. Each fish was then rapidly transferred back to 14°C for recovery. The use of CT_max is a standard approach to assess upper thermal tolerances in ectotherms and is particularly useful for comparisons across species (Beitinger et al., 2000).

To assess the transcriptome-wide effects of exposure to an environmentally relevant temperature, we exposed a different subset of ∼50 dph longfin smelt and delta smelt to 20°C. Water temperature for eight fish at 14°C was raised at a rate of 0.3°C min⁻¹ until it reached 20°C. Once the target temperature was reached, the fish were rapidly transferred back to 14°C and allowed to recover from exposure for 60 min (n=8 in each treatment). This acute exposure protocol allows for the detection of a strong cellular stress response in delta smelt at temperatures that reflect their upper thermal distributions in the wild (Komoroske et al., 2014, 2015). The Delta also undergoes diurnal temperature fluctuations that result in resident fishes experiencing rapidly changing temperatures daily (Komoroske et al., 2015). A 14°C acclimation temperature treatment group that was treated identically as the 20°C group was used as a control. After 60 min of recovery, fish were killed in buffered tricaine methanesulfonate, frozen in liquid nitrogen and stored at −80°C. Handling, sample sizes and temperature exposures were identical for both species.

Metabolic rates (Ṁ_O2) for delta smelt and longfin smelt were determined by measuring oxygen consumption at 14°C and 20°C in a four-chamber, flow-through respirometry system allowing for multiple closed-chamber Ṁ_O2 measurements for each individual. Oxygen consumption was measured by recording changes in air saturation with a fiber-optic oxygen meter (Witrox 4; Loligo Systems, Tjele, Denmark) and individual oxygen sensor spots (accuracy of ±0.4% at 20.9% O₂) fixed within each 7 ml chamber (Loligo Systems). Chambers were held in a temperature-controlled circulating freshwater bath that was continuously mixed and aerated. Trials consisted of three chambers containing a single fish and a control chamber to account for biological activity in the water. Fish were individually transferred to a chamber and held for a 20 min habituation period at 14°C post-transfer. For the 14°C treatment, each chamber was flushed with aerated water following the habituation period until the chamber returned to ∼100% air saturation; the chamber was then closed and oxygen consumption was recorded until at least a 5–10% decrease in air saturation was reached. The cycle of flushing and sealing the chamber for oxygen consumption measurements was repeated two to five times using a peristaltic pump (iPump i150; Loligo Systems). Measurements for the 20°C treatment were similar to the 14°C treatment; however, following the habituation period, the water bath temperature was increased to 20°C at a rate of 0.3°C min⁻¹ before oxygen consumption was measured. No chamber went below 75% air saturation and no supersaturation (>100% air saturation) occurred while warming the water. Following each trial, the wet mass of the fish was measured. Although fish may need up to 24–48 h after transfer and handling stress to allow Ṁ_O2 to return to a resting state (Clark et al., 2013), this was not possible in the present study because of the sensitivity to isolation and starvation for these species during early development in captivity (Lindberg et al., 2013). Additionally, these fish are sensitive to light exposure; consequently, the Ṁ_O2 treatments were conducted in the dark, preventing observations of changes in activity. Therefore, although true resting Ṁ_O2 values are difficult to obtain with these species, the goal was to conduct a comparative study of the effects of water temperature on two sensitive species that were treated identically. Because no Ṁ_O2 measurements exist for this early developmental stage, these data are the first estimates of Ṁ_O2 for larval delta smelt and longfin smelt.

We quantified Ṁ_O2 using an approach modified from Steffensen et al. (1994), described in Clark et al. (2013). Briefly, the first 20 min was excluded to eliminate noise in the measurements associated with handling stress and flushing/sealing the chamber. The remaining closed-chamber Ṁ_O2 measurements were calculated using two-point regressions (delta smelt, 14°C n=14, 20°C n=14; longfin smelt, 14°C n=13, 20°C n=9). A frequency distribution of regression values was generated, yielding a bimodal distribution in each temperature treatment. One cluster of regressions included the control chambers used for background activity, while a second distribution was for actively respiring fish. Background activity in the control chambers were averaged for each temperature treatment and subtracted from the respective Ṁ_O2 measurements. No differences in Ṁ_O2 measurements were detected between repeated measurements (ANOVA) for each individual; therefore, to facilitate comparisons between treatments and species, all Ṁ_O2 measurements used for statistical analyses were taken at a similar time point post-transfer into the chamber [mean time of measurement was 48.8 min (±3.2 min s.e.m.) post-transfer]. The final Ṁ_O2 was corrected for fish mass and is presented as µmol O₂ h⁻¹ g⁻¹.

Whole fish were homogenized using a Qiagen TissueLyser in Buffer RLT Plus (RNeasy Plus Mini Kit); 350 μl of the homogenate was used for RNA extraction on a QIAcube following the manufacturer's protocols (Qiagen, Valencia, CA, USA). Five samples from the 14°C and 20°C treatments from each species (n=20 fish total) were selected for sequencing. The RNA quality was assessed using a Bioanalyzer (Agilent, Santa Clara, CA, USA; RNA integrity numbers=9.2–10 for all samples).

Sequencing was performed at the Vincent J. Coates Genomic Sequencing Laboratory at the University of California, Berkeley (supported by National Institutes of Health S10 Instrumentation grants S10RR029668 and S10RR027303). Total RNA was used to prepare separate Poly-A isolated, strand-specific cDNA libraries using Directional PrepX RNA Library Prep kits (Wafergen Biosystems, Fremont, CA, USA) prior to sequencing. Individual samples were barcoded with unique adapters supplied by the manufacturer. Because species-specific reference transcriptomes were required for this study, we conducted two separate differential expression experiments. Therefore, libraries for each species were pooled and sequenced on an Illumina HiSeq 2500 in rapid run mode over two lanes (i.e. two lanes for each species). An additional seven longfin smelt (collected from the wild or from previous years) were sequenced over two lanes to improve transcriptome coverage for the assembly. Because there was no previous transcriptome information for longfin smelt, sequencing was conducted to produce 100 base pair paired-end reads, generating a mean (±s.e.m.) of 33.1±1.3 million raw reads per individual. A delta smelt transcriptome had previously been assembled from reads generated on the 454 sequencing platform (Jeffries et al., 2015). Therefore, 100 base pair single-end reads were generated for delta smelt, producing a mean (±s.e.m.) of 36.6±2.6 million raw reads per individual.

Raw, de-multiplexed reads were filtered and trimmed for adapter contamination and low quality sequences (using the programs Scythe and Sickle, respectively; available at https://github.com/ucdavis-bioinformatics/) to generate means (±s.e.m.) of 30.3±1.2 million and 33.5±2.4 million trimmed reads for longfin smelt and delta smelt, respectively. The longfin smelt reads were normalized and assembled using Trinity v.r20140413p1 (Haas et al., 2013) to produce a de novo transcriptome assembly with 125,843 sequences (contigs). The raw assembly was filtered to remove rRNA, low or unexpressed transcripts, and transcripts with no annotation [e.g. a blast hit, sequence description or gene ontology (GO) term] to a total of 69,949 transcript contigs, representing 42,954 genes.

Based on an initial alignment of the Illumina reads for the delta smelt, the existing reference transcriptome was not comprehensive enough to include all genes expressed in the present experiment. Therefore, all trimmed reads were first aligned to the 454 isotig set from the previous transcriptome assembly (generated using Newbler software; Roche, Pleasanton, CA, USA; details in Jeffries et al., 2015) with Burrows–Wheeler Aligner (BWA) short read aligner v.0.6.2 (Li and Durbin, 2010). The reads that did not align were normalized and assembled using Trinity to generate 67,556 contigs. These two isotig/contig sets were concatenated (101,051 sequences). This combined assembly was filtered the same way as with the longfin smelt assembly to a total of 60,737 transcript contigs representing 36,788 genes.

The transcript contigs for each assembly were annotated by first blasting each contig to a 521,381 sequence subset of the National Center for Biotechnology Information (NCBI) nr protein database representing Actinopterygii (ray finned) fishes. Supplemental blasts were performed against the NCBI swissprot and nr databases to enhance annotation. Up to 20 ‘top hits’ (e-value cutoff of 1E–05) from the blast results were input to Blast2GO for functional annotation (Conesa et al., 2005; Götz et al., 2008), which was also used to assign GO terms to each transcript.

The reads were then aligned to the appropriate reference transcriptome using BWA. For each species (n=10 species⁻¹), a table of raw counts by transcript (from the BWA alignments to the filtered assembly) was generated using sam2counts.py (https://github.com/ucdavis-bioinformatics). The read counts for each transcript contig were summed to gene level to generate a table of raw counts by gene for the differential expression analysis, which was conducted using edgeR (Robinson et al., 2010). Genes were considered differentially expressed at a Benjamini–Hochberg corrected false discovery rate (FDR) of <0.05. Because both species were sequenced separately and reads were mapped to different assembled transcriptomes, the gene-by-gene response patterns are not directly comparable and therefore were analyzed separately for each species. To facilitate comparisons of the responses between these species, we grouped the differentially expressed transcripts together based on the GO terms (level 2). The GO terms represent the predicted gene function(s) assigned to each contig. We used a functional analysis approach in Blast2GO for comparing the broad patterns in the gene expression profiles with a larger significant gene list (FDR<0.1). Functional GO categories were considered enriched in the significant gene list at P<0.01 (Fisher's exact test). Complete lists of differentially expressed transcripts at FDR<0.1 for both species are provided in the supplementary material (Tables S1, S2). The raw sequence data are available through the NCBI Sequence Read Archive with the accession number SRP064394.

Every month from December 2013 to July 2015 except one has been below −4 on the Palmer Drought Severity Index, indicating extreme drought (Fig. 1A) and providing evidence that this was the driest period in California from 1995 to 2015. During the drought years of 2012–2015, mean monthly water temperatures in the freshwater regions of the Delta from April to July were on average higher than between 1995 and 2011 (two-factor ANOVA, Bonferroni corrected P<0.01), demonstrating the effects of drought on surface water temperatures (Fig. 1B). Wild delta smelt on average (18.9±2.35°C, mean±s.d.) were found at warmer temperatures than wild longfin smelt (16.4±2.35°C) between March and July (1995–2015; Fig. 2). This corresponded with the overall upper thermal tolerances of the two species as the CT_max for delta smelt was 27.6°C (±0.28°C s.e.m.), which was significantly higher than the longfin smelt CT_max of 24.8°C (±0.38°C s.e.m.; Mann–Whitney U-test, P<0.001; Fig. 2).

Oxygen consumption was higher overall in 50 dph delta smelt than in longfin smelt (two-factor ANOVA, P<1.5×10^–7; Fig. 3). The mean (±s.e.m.) Ṁ_O2 for delta smelt at 20°C (65.5±5.1 µmol O₂ h⁻¹ g⁻¹) was significantly higher (P<0.05) than that at 14°C (47.9±5.0 µmol O₂ h⁻¹ g⁻¹). Conversely, there was no statistically significant difference in Ṁ_O2 for longfin smelt at 20°C and 14°C (P=0.77). Mean (±s.e.m.) Ṁ_O2 at 20°C was 21.9±3.0 µmol O₂ h⁻¹ g⁻¹ and that at 14°C was 29.0±4.5 µmol O₂ h⁻¹ g⁻¹ for longfin smelt.

The transcriptomic profiles for both species indicated that acute exposure to 20°C resulted in the differential expression of genes involved in growth and development, in addition to genes involved in ion regulation. Additionally, there were species-specific transcriptomic responses to elevated temperature. In longfin smelt, there was differential expression of many genes associated with thermal stress, which included the upregulation of heat shock proteins and protein chaperones that are generally characteristic of a heat shock response in fishes (Basu et al., 2002). Conversely, delta smelt responses involved the differential expression of genes associated with protein synthesis and metabolic processes. The differences in the suite of genes that were differentially expressed likely reflect the relative magnitude of the temperature stressor between the two species.

There were 229 transcripts differentially expressed in the longfin smelt due to the temperature treatments at FDR<0.1 used for the functional analysis. Several GO terms associated with growth and development were enriched in the significant gene list (Fig. 4A). These included categories directly linked with development (e.g. cartilage development involved in endochondral bone morphogenesis, branching morphogenesis of a nerve), but also those associated with collagen stabilization (collagen catabolic process) and the extracellular matrix (e.g. extracellular matrix disassembly, extracellular matrix structural constituent) that are precursors to bone growth and skeletal development (Vieira et al., 2013). Several GO terms that could be associated with responding to external stimuli (growth factor binding, response to activity, signal transduction in absence of ligand, regulation of extrinsic apoptotic signaling pathway) were also enriched in the significant gene list. Lastly, the GO term ion channel activity was significantly enriched, indicating a potential disruption of cellular ion homeostasis.

At FDR<0.05, 92 transcripts were upregulated and 60 transcripts were downregulated (Table S1) in longfin smelt exposed to 20°C compared with the 14°C group. Many of the upregulated genes are involved in a response to temperature and a general stress response (Fig. 5A). Consistent with this pattern was the upregulation of protein chaperones (heat shock protein beta-11-like, heat shock protein 27 beta and dnaJ homolog subfamily C member 21-like) and transcripts involved in aerobic metabolism (citrate synthase, mitochondrial-like and carbonic anhydrase 14-like) and response to DNA damage (DNA-dependent protein kinase catalytic subunit).

There was also differential regulation of genes involved in cell signaling pathways, consistent with the functional analysis. Of the signaling genes impacted by the temperature treatments, several genes involved in the MAPK/ERK signaling pathway were upregulated, including dual specificity protein phosphatase 7, which is involved in regulating MAPK activity in response to changes in the cellular environment. There was also a downregulation of kinase suppressor of RAS 2-like, which inhibits the MAPK (Channavajhala et al., 2003) regulation of the ERK, JNK and NF-kappa-B signaling pathways. The NF-kappa-B signaling pathway was impacted as there was an upregulation of mitogen-activated protein kinase kinase kinase 14-like and a downregulation of inhibitor of nuclear factor kappa-B kinase subunit alpha-like and inhibitor of Bruton tyrosine kinase-like. Lastly, transforming growth factor-beta receptor type I b, which is involved in many cell signaling pathways that are involved in cell proliferation and growth, was upregulated at 20°C.

Growth and development genes were differentially expressed at 20°C and many genes were associated with the extracellular matrix, biological adhesion, cell junction, growth and structural activity (Fig. 5), which is consistent with the functional analysis. There were several collagen genes that were upregulated at 20°C [e.g. collagen alpha-1(II) chain-like isoform 1, collagen alpha-1(IX) chain-like and collagen alpha-1(XII) chain-like]. Additionally, genes involved in bone growth were affected by the temperature treatments: periostin precursor and tensin-3-like were significantly upregulated whereas fibroblast growth factor receptor 2-like isoform X8 was downregulated at 20°C.

Numerous genes associated with calcium signaling and muscle activity (i.e. locomotion; Fig. 5A) were impacted in the longfin smelt at 20°C. Associated with muscle activity is the regulation of cellular calcium and potassium levels; however, these processes can also be affected by osmoregulatory stress (Brennan et al., 2015). Longfin smelt showed differential expression of several genes involved in ion channels and ion homeostasis due to the temperature treatment (e.g. aquaporin-4 and transient receptor potential cation channel subfamily M member 4-like were significantly upregulated at 20°C, whereas voltage-gated potassium channel subunit beta-3-like, Kv channel-interacting protein 4-like isoform X1 and voltage-dependent L-type calcium channel subunit alpha-1F-like were downregulated), a pattern supported by the functional analysis results.

There were 243 transcripts differentially expressed at FDR<0.1 in the delta smelt used for the functional analysis (Fig. 4B). In contrast to the patterns detected in longfin smelt, many of the GO terms enriched in the delta smelt significant gene list were associated with protein synthesis (e.g. RNA localization, transcription initiation from RNA polymerase II promoter) and metabolic processes (glucose metabolic process, aspartate family amino acid metabolic process). There was also an enrichment of GO terms associated with development (e.g. cell fate commitment, establishment of tissue polarity).

There were 95 upregulated and 58 downregulated transcripts at FDR<0.05 in delta smelt exposed to 20°C compared with the 14°C group (Table S2). Many of the genes differentially expressed were consistent with a response to changes in water temperature; however, no protein chaperones or heat shock proteins were upregulated at 20°C. Several of the differentially expressed genes are involved in transcription regulation and protein synthesis, which included the upregulation of inducible transcription factors (immediate early response gene 2 protein, fos-related antigen 2-like, transcription factor jun-B-like, transcriptional regulator Myc-2 and proto-oncogene c-fos). There was also an upregulation of genes involved in cell proliferation and DNA repair (serine/threonine-protein kinase PLK2-like isoform 2 and DNA repair protein RAD51 homolog A). Many genes involved in general metabolism (Fig. 5A) were impacted, including genes involved in fatty acid metabolism (acetyl-CoA carboxylase 1-like, hydroxymethylglutaryl-CoA synthase, cytoplasmic-like isoform X1 and 3-hydroxyacyl-CoA dehydratase-like) and polyamine synthesis (ornithine decarboxylase-like), involved in regulating cellular processes such as growth and proliferation (Larqué et al., 2007). However, creatine kinase-2 and alpha amylase, both involved in metabolic processes (Connon et al., 2011), were downregulated.

Similar to the longfin smelt patterns, many genes associated with growth and development were differentially expressed in the delta smelt at 20°C; however, fewer genes were associated with the extracellular matrix (Fig. 5B). Genes involved in regulating the Wnt signaling pathway (low-density lipoprotein receptor-related protein 6-like and protein naked cuticle homolog 1-like), a pathway that is involved in regulating chondrocytes and osteoblasts (Vieira et al., 2013), and in muscle function and development (Cisternas et al., 2014), were differentially expressed at 20°C. Genes involved in ion regulation were also differentially expressed in the 20°C treatment (probable cation-transporting ATPase 13A3-like was upregulated, whereas sodium/potassium/calcium exchanger 6-like and sodium leak channel non-selective protein precursor were downregulated).

Similar to the pattern found in longfin smelt, there was an activation of the NF-kappa-B signaling pathway with an upregulation of TNF receptor-associated factor 5, which mediates the NF-kappa-B signaling pathway, and the downregulation of SET domain-containing protein 6 and TNFAIP3-interacting protein 1-like, which are associated with repression of the NF-kappa-B transcription factor activity. Several genes that are associated with immune response (growth factor receptor-bound protein 10-like and toll-interacting protein-like) and MAPK (dual specificity protein phosphatase 5-like and transmembrane protein 184B-like) signaling were also impacted.

We conducted the first comparative study on the cellular and physiological responses to elevated temperature of two estuarine fishes of conservation concern. Our study demonstrates the importance of integrating transcriptome-wide responses with whole-organism responses to identify the effects of high water temperature on sensitive species. The CT_max estimates are close to the uppermost temperatures where these fish are found in the wild (0.8°C and 0.6°C higher for longfin smelt and delta smelt, respectively), suggesting that these fish may already be living close to their thermal limits. Given the limited distribution of the delta smelt and the lack of evidence indicating exploitation of more saline habitats (Bennett, 2005), it is unlikely that the delta smelt will be able to adjust their range in response to future warming or periods of extreme drought. In contrast, longfin smelt exhibit a more anadromous and euryhaline life history, suggesting the potential to move to suitable habitats by shifting their distribution between the estuary and the open ocean or to other estuaries if necessary.

Our data suggest that chronic exposures to temperatures ≥20°C may be unsustainable for juvenile longfin smelt. The CT_max for longfin smelt relative to the delta smelt suggests a lower thermal tolerance for this species, which was consistent with the transcriptome-wide responses at 20°C. Additionally, we did not detect a change in oxygen consumption in longfin smelt over a 6°C increase in water temperature despite an upregulation of genes involved in aerobic metabolism. A mismatch between whole-organism and cellular-level aerobic metabolism could arise for several reasons. It may suggest that the longfin smelt did not have sufficient steady-state proteins to cope with the increased metabolic costs and the upregulation of citrate synthase during recovery represents a cellular shift towards responding to increased oxygen demand. However, fishes may have a limited metabolic scope during larval stages, which can impact their ability to tolerate multiple stressors (Killen et al., 2007). This potentially supports our observation as the cumulative effects of confinement, temperature and handling decreased our ability to detect a temperature-dependent increase in oxygen consumption in longfin smelt. It is important to note that the sampling times were different for the fish used for the respirometry trials compared with those used for the transcriptome-wide analyses, limiting our ability to determine the mechanistic cause of the different response patterns with certainty. The longfin smelt patterns are in contrast with those in delta smelt, where we detected increased oxygen consumption to match a presumed increase in demand at the cellular level. The increase in the metabolic rate of delta smelt occurred at a temperature that resulted in elevated expression of genes involved in overall metabolism, which is consistent with a greater metabolic demand at high temperatures in ectotherms.

Many of the genes affected by elevated temperature for both species were associated with muscle function, growth and skeletal development; responses also observed in muscle and skeleton transcriptome sequencing studies on other fishes (Long et al., 2012; Palstra et al., 2013; Vieira et al., 2013). High water temperature can alter expression of genes involved in extracellular matrix and cytoskeleton restructuring, and collagen stabilization in fishes (Jeffries et al., 2012a; Liu et al., 2013; Tomalty et al., 2015), processes that are critical during skeletal growth (Vieira et al., 2013) and were impacted in the present study. Because the larval stages of development in fishes are plastic, temperature impacts during early life stages may lead to long-term effects on growth (Johnston, 2006). It is important to note that because of the small size of the fish, whole-body homogenates were used for the transcriptome-wide analyses to ensure sufficient amounts of RNA for sequencing. Consequently, we increased the likelihood of detecting genes associated with growth and development, and potentially limited our detection of tissue-specific expression patterns. However, this was the first transcriptome sequencing for longfin smelt and it was to build on the transcriptome information available for delta smelt; therefore, we determined beforehand that it was best to generate the most comprehensive and representative set of sequences for de novo transcriptome assembly.

We detected the differential expression of genes involved in ion regulation in longfin smelt and delta smelt, suggesting that osmoregulatory function might be impacted at elevated temperatures. Previous work has shown alterations in gene expression, enzyme activity and blood plasma ion levels associated with osmoregulation in adult Pacific salmon exposed to elevated temperatures (Crossin et al., 2008; Jeffries et al., 2012b, 2014). This could be a critical indirect consequence of exposure to elevated temperature in anadromous or estuarine fishes. Impacts on osmoregulatory ability or preparedness for saltwater conditions may delay migration and impact estuary survival, either directly or indirectly through increased susceptibility to predation and disease (McCormick et al., 2009; Halfyard et al., 2013). In the Delta, salinity levels can vary dramatically in response to tidal fluctuations and rivers discharging into the estuary (Cloern et al., 2011). During periods of drought, reduced river flows may lead to increased saltwater influx into the Delta. Because these periods coincide with high water temperatures, it is possible that alterations in osmoregulatory ability associated with temperature increases may impact longfin smelt and delta smelt during drought years.

Consistent with a lower thermal tolerance, longfin smelt exhibited a pronounced cellular stress response after exposure to 20°C that was not observed in delta smelt. Cellular responses (and their onset) are more sensitive to changes in temperature than whole-organism responses, such as upper thermal tolerance limits (Chadwick et al., 2015), and are indicative of temperatures that begin to adversely affect an organism (Kültz, 2005; Kassahn et al., 2009). Longfin smelt exhibited a cellular stress response at a temperature below those at which they can be periodically found in the wild, suggesting that these fish may routinely experience sublethal temperature stress in the Delta. There was an upregulation of heat shock proteins, protein chaperones, and genes associated with DNA damage and aerobic metabolism, which are all characteristic of an acute temperature stress response in ectotherms (Kassahn et al., 2009). The upregulation of heat shock protein beta-11-like and heat shock protein 27 beta in longfin smelt was consistent with an acute thermal stress response in juvenile hybrid catfish (Liu et al., 2013). Many of the genes affected by the temperature treatment are also involved in important developmental processes (Long et al., 2012). Because the expression patterns are being compared with those of a control group at the same developmental stage, we can attribute these changes to exposure to 20°C and not to developmental differences. However, the present results are likely confounded by an interaction between a temperature response and normal developmental processes.

The gene expression patterns also suggest the activation of the MAPK/ERK and NF-kappa-B signaling pathways. The ERK signaling pathway has a role in cell division, migration and survival and is often stimulated by extracellular growth factors (Cowan and Storey, 2003), consistent with the differential expression of genes associated with growth factor binding in longfin smelt. In the fine flounder (Paralichthys adspersus), extracellular growth factors and MAP/ERK signaling were shown to be crucial for skeletal muscle and overall growth (Fuentes et al., 2011). We also detected an upregulation of mitogen-activated protein kinase kinase kinase 14-like, which is involved in regulating the NF-kappa-B signaling pathway (Secombes, 2008), consistent with patterns observed in Pacific salmon (Jeffries et al., 2014). The activation of the NF-kappa-B signaling pathway observed in longfin smelt was also characteristic of an acute temperature stress response in rainbow trout (Oncorhynchus mykiss; Lewis et al., 2010). The NF-kappa-B signaling pathway has also been associated with an immune and inflammation response in fishes (Mu et al., 2010). Given that many genes associated with an inflammation response are also related to restructuring the extracellular matrix, the activation of the NF-kappa-B signaling pathway in the present study may also be linked with the differential expression of genes involved in extracellular matrix processes and skeletal development.

The 20°C temperature treatment in the present study was well below the upper thermal limits of the delta smelt (∼8°C below, Komoroske et al., 2014; 7.6°C below in the present study) and within the range of temperatures this species experiences in the wild. Accordingly, exposure to 20°C resulted in the altered regulation of a different suite of genes than longfin smelt, with many of the differentially expressed genes in delta smelt involved in metabolic processes and protein synthesis. This likely reflects differences in the relative magnitude of the thermal stressor between the two species as 20°C is further from the upper thermal limits of delta smelt. We detected the upregulation of several inducible transcription factors that belong to a group of genes known as immediate early genes. These inducible transcription factors are relatively conserved across taxa (Kültz, 2005) and are associated with regulating genes that promote cell survival after exposure to changes in environmental conditions (Kassahn et al., 2009). Interestingly, fos-related antigen 2-like has been associated with osteoblast function and bone and extracellular matrix development (Bozec et al., 2010); therefore, it may be linked to the effects of temperature on growth and skeletal development in delta smelt. There were many genes associated with muscle function that were differentially expressed, potentially because of the larger size of the delta smelt relative to the longfin smelt at 50 dph. Some of the functional categories that were enriched in delta smelt were similar to those in larval zebrafish (Danio rerio) exposed to elevated temperatures (Long et al., 2012). This suggests that some of the patterns observed in delta smelt may be part of a general response to elevated temperature in larval fishes.

The upregulation of inducible transcription factors has been suggested to be involved in an organism's cellular response as temperatures approach upper thermal limits (e.g. Kassahn et al., 2009). However, we detected an upregulation of inducible transcription factors at temperatures 4–6°C below those that cause an extreme heat shock response in delta smelt (Komoroske et al., 2015). Interestingly, during periods of extreme temperature stress, a different inducible transcription factor is upregulated in delta smelt (i.e. AP-1; Komoroske et al., 2015) than those detected in the present study at a more moderate temperature exposure. Inducible transcription factors were upregulated during relatively routine increases in water temperature associated with tidal fluctuations that were well below the upper thermal limits of the intertidal California mussel (Mytilus californianus; Connor and Gracey, 2011). In addition to activating the MAPK and NF-kappa-B signaling pathways that maintain normal cell function (Cowan and Storey, 2003; Kassahn et al., 2009), the upregulation of immediate early inducible transcription factors may promote an adaptive response to environmental stress (Cohen, 1997). Having rapidly inducible cellular responses to environmental changes may be consistent with the ability to cope with fluctuating local environmental conditions, common in dynamic estuary ecosystems.

We used an integrative approach to assess the impact of elevated water temperature on two fishes of immediate conservation concern in a large estuary ecosystem. The environmentally relevant exposure temperature of 20°C is more commonly reached and occurs earlier in the year during periods of extreme drought. Based on their physiological responses, longfin smelt may be more susceptible to increases in temperature, as this species appears to be near its upper thermal limits with relatively little room to tolerate persistent drought or future climate warming in California. Examination at multiple levels of biological organization (cellular to whole-organism) highlighted the potential vulnerability of longfin smelt relative to delta smelt. Longfin smelt aggregate in deeper water during peak seasonal temperatures, potentially to avoid extreme surface temperatures (Rosenfield and Baxter, 2007). Therefore, longfin smelt may respond behaviorally by seeking out suitable temperatures. Alternatively, longfin smelt may need to adjust the phenology of their spawning or migration to the cooler saltwater environment; otherwise, periods of high water temperature could be detrimental to segments of this population as regions of the Delta will no longer provide suitable habitat. With the delta smelt at a legitimate risk of extinction, longfin smelt could become a key indicator species of ecosystem health in the Delta, the epicenter of the debates surrounding water use in California. Understanding the physiological responses to environmental stressors of threatened and endangered species is crucial for conservation efforts in coastal systems and to effectively manage important estuary ecosystems.

We are grateful to Galen Tigan, Luke Ellison, Joan Lindberg and Tien-Chieh Hung at the University of California, Davis Fish Conservation and Culture Laboratory for spawning and rearing the fish prior to the exposures and for helpful discussions. We thank Randy Baxter and Jennifer Messineo at the California Department of Fish and Wildlife and Matthew Dekar at the United States Fish and Wildlife Service for assistance in sample collections. We also thank Jennifer Truong for assistance with processing samples, Blythe Durbin-Johnson for advice on transcriptomics analyses, Dennis Cocherell for technical advice, and Minyong Chung and the Vincent J. Coates Genomics Sequencing Laboratory at the University of California, Berkeley. We also thank an anonymous reviewer for helpful comments on an earlier version of this manuscript.