Comparative physiological, biochemical and molecular thermal stress response profiles for two unionid freshwater mussel species

Over the past decade, global mean surface air temperatures have increased by 0.81°C and are predicted to increase by an additional 1.4–3.1°C by the year 2100 (IPCC, 2014), with the southeastern USA falling within the middle of this range (Burkett et al., 2001; van Vliet et al., 2013). As water temperatures tend to parallel air temperatures, freshwaters have also experienced a warming trend. In the southeastern USA, summer freshwater temperatures can average 35°C and exceed 40°C during extreme heat waves, which is near or above the upper thermal tolerance limits of many freshwater species (Archambault et al., 2014; Pandolfo et al., 2010; Ganser et al., 2015; also reviewed by Ficke et al., 2007). Such significant temperature changes have the potential to disrupt ecosystem stability, resulting in species-specific differences in susceptibility, with factors such as physiological tolerance range, potential for adaptation to abiotic stressors, and life history characteristics (Hofmann and Somero, 1995; Somero, 2002, 2012) all playing important roles in determining a ‘winners’ versus ‘losers’ status (Somero, 2010).

A substantial portion of all freshwater mussel species, including the considerable diversity found within the southeastern USA, are in global decline (Strayer, 2006). These declines most likely result from a combination of physical habitat modification, point source pollution, sedimentation, invasive species, and warming temperatures due to climate change (as reviewed by Hastie et al., 2003). Warming trends have already caused alterations in thermal structure and water quality (Cruise et al., 2000; reviewed by Caissie, 2006), resulting in downstream effects on primary producers, migratory fishes and macroinvertebrates (Davis et al., 2013; Harper and Peckarsky, 2006; MacFadyen et al., 2004; O'Reilly et al., 2003; Winder et al., 2009). Ecosystem-level functional alteration may be especially pronounced within the state of Alabama, a renowned biodiversity hotspot for freshwater fauna, particularly fishes, crayfish and mollusks (Lydeard and Mayden, 1995). Mussels provide fundamental services such as nutrient cycling (Christian et al., 2008; Jost and Helmuth, 2007), habitat modification (Norkko et al., 2006) and water filtration (Budd et al., 2001). This functional diversity thus supports increased biodiversity within aquatic systems (Donadi et al., 2013; Vaughn et al., 2008). Globally, there are approximately 670 species of freshwater mussel, of which 300 can be found in North America (Cummings and Graf, 2010). Around 182 freshwater mussel species have been identified in the state of Alabama, more than in any other state in the USA (Williams et al., 1993; ADCNR, 2015). However, 26 are considered extinct (Haag, 2012) and 57 are currently listed as of high conservation concern according to Alabama's Wildlife Action Plan (ADCNR, 2015). Such biodiversity loss is particularly disastrous for Alabama's aquatic ecosystems because freshwater mussels are keystone species and future extirpations will lead to ecosystem health declines (Spooner and Vaughn, 2006; Vaughn et al., 2008).

A significant amount of climate change-related research on bivalve mollusks has focused on physiological variation in response to thermal stress within and between marine bivalve species distributed along latitudinal gradients (Somero, 2005; Compton et al., 2007; Morley et al., 2009a,b; Jones et al., 2009; Chauvaud et al., 2012). For freshwater mussel species, there are three general types of distribution pattern within the USA; species that are widely distributed across latitudinal biogeographic gradients are contrasted by those that are restricted to either northern or southern geographic distributions as a result of warm or cool temperature limitations, respectively (Haag, 2012). Little focus has been placed on understanding the inherent differences in physiological robustness that allow one species to widely expand their geographic range while another congener species remains restricted to a smaller region. Of particular interest are questions regarding mechanisms that determine ‘winner’ versus ‘loser’ status when the geographic range of a narrowly distributed species exists within that of a widespread congener species. Under the threat of climate change, are widespread species more likely to be ‘winners’ and narrowly distributed congener species more likely to be ‘losers’? Such research can provide informative data required for successful conservation and management of freshwater mussels under the continuing threat of global climate change.

Previous research into how freshwater mussels might respond to future global warming has focused on acute heat shock and basic mechanistic experimental designs (Strayer and Dudgeon, 2010; Strayer, 2006). However, integrative experimental approaches have also successfully yielded results more applicable to conservation and management (Dégremont et al., 2007; Clark et al., 2001; Samain et al., 2007). Furthermore, correlating gene expression data with physiological and traditional life history parameters in response to chronic sub-lethal heat stress has become more prevalent. Such integrative approaches have been embraced across multiple disciplines and research interests, including freshwater (Doucet-Beaupré et al., 2010; Wang et al., 2012) and marine bivalves (Lang et al., 2009; Bettencourt et al., 2010; Clark et al., 2010; Craft et al., 2010; Fleury et al., 2010).

To implement an integrative approach for assessing thermal tolerance in freshwater mussels, a series of experiments were performed to assess the molecular (candidate gene expression), metabolic (glycogen and triglyceride content) and physiological (mortality, condition index, growth rates) responses of two unionid freshwater mussel species, Villosa lienosa and Villosa nebulosa, under environmentally relevant acute heat shock and chronic thermal warming scenarios predicted for 2100 in the southeastern USA. These two species were chosen based on their close evolutionary history, as well as differences in conservation status and current geographic range. Specifically, V. lienosa is widespread throughout the southeastern USA and very stable within its range. In contrast, V. nebulosa is of high conservation concern and geographically restricted to the Coosa, Cahaba and Black Warrior rivers in the Mobile River Basin above the fall line within the states of Alabama, Georgia and Tennessee, within the middle of the geographic range of the broadly distributed V. lienosa (Fig. 1). Based on differences in geographic range, we hypothesized that V. lienosa was the more physiologically robust species, and this difference in robustness would be represented by a greater thermal tolerance. By measuring mortality, gene expression, condition index, growth rates, and glycogen and triglyceride content, we investigated differential response mechanisms to a chronic warming simulation that could be linked to each species' unique thermal tolerance range. Along with presenting conclusive evidence that V. lienosa is more thermally tolerant than V. nebulosa, we provide further evidence that future global warming will impact current freshwater mussel populations in a species-specific manner. Our results also support the argument that long-term experiments in natural field environments using realistic thermal regimes can provide informative data required for successful conservation and management of freshwater mussels.

Three adult V. lienosa (Conrad 1834) were collected in May 2011 from wild populations within Sturdivant Creek, Monroe County, AL, USA. Seven V. nebulosa (Conrad 1834) were collected in April 2011 within the South Fork of Terrapin Creek and Shoal Creek, Cleburne County, AL, USA. Glochidia were harvested from charged females and cultured using standard practices at the Alabama Aquatic Biodiversity Center (AABC), a department of the Alabama Department of Conservation and Natural Resources recovery and research program based in Marion, AL, USA. Transformed progeny were cultured for 12–16 months. Acute temperature challenges were performed with ∼14 month old post-transformed cultured mussels and chronic temperature challenges were performed with ∼16 month old post-transformed cultured mussels.

Adult V. lienosa and V. nebulosa were acclimated to laboratory temperatures of 22°C for 17 days in 10 gallon (∼38 l) tanks with daily feedings of Shellfish Diet 1800^® (Reed Mariculture) and water changes every 3 days with AABC site water taken from aquaculture ponds. Nine mussels were sampled at 0 h for basal expression, and the remaining mussels were placed in 4 l Nalgene beakers containing 3500 ml of AABC site water that was acclimated to either ambient (22°C) or heat challenge (30°C) temperature. Each Nalgene beaker contained 12 mussels and three beakers were used for each species and temperature treatment (12 beakers total). Three mussels were sampled from each beaker (total of nine mussels per treatment) at 3, 6 and 24 h during the heat shock and after 48 h of recovery at 22°C. Whole tissue was immediately dissected, flash frozen in liquid nitrogen and stored at −80°C. Total RNA was extracted from 50 mg frozen tissue using TRIzol^® Reagent (Invitrogen™) according to the manufacturer's instructions, with the inclusion of the additional high-salt precipitation. Three biological replicates for each species and treatment were generated by pooling equivalent amounts of total RNA from the three individual mussels sampled from each beaker at each time point.

The iSCRIPT™ cDNA kit (Bio-Rad, Hercules, CA, USA) was used to reverse transcribe 1 μg of pooled total RNA to cDNA for gene expression analysis. Three genes were selected from a previously constructed 454 Ion Torrent transcriptomic database (Perkins and Jenny, 2012) based on functional annotation to classic heat shock proteins (HSP90, HSP70-1, HSP70-2). Primers for reverse transcription real-time quantitative PCR (RT-qPCR) are listed in Table 1. qPCR was performed using the PerfeCTa^® SYBR^® Green Supermix for iQ™ (Quanta Biosciences, Gaithersburg, MD, USA) in a MyiQ2 Two-Color Real-Time PCR Detection system (Bio-Rad). qPCR conditions used were: 95°C for 3 min, 94°C for 15 s and 57°C for 30 s (40 cycles). At the end of each PCR run, a melt curve analysis was performed to ensure that only a single product was amplified. Three technical replicates were performed for each biological replicate. For quantification of gene expression, a standard curve was generated by serially diluting plasmids containing a fragment of the target transcript. All transcripts were normalized to the ribosomal 18S gene. The ratio of induction was determined for each replicate by dividing the heat shock transcript number in order to estimate fold changes.

Statistical analysis of gene expression data was performed using SAS software (version 9.3, Cary, NC, USA) and all graphs were produced using the Prism 5 software (GraphPad Software Inc., San Diego, CA, USA). Acute heat shock gene expression profiles were analyzed via two-way ANOVA (α=0.05). Basal gene expression profiles were analyzed via unpaired t-tests (α=0.05). Values for overall and basal expression are presented as means±s.d., and induction ratios are presented as means±s.e.m. due to repeated comparisons to the average transcript number for every time point.

To simulate future warming in a natural and variable environment, two established field ponds located at the AABC were maintained at ambient temperature (un-manipulated standard thermal aquaculture conditions) or approximately 3°C above the ambient pond temperature via reduction of the constant flow of artesian well water into the pond, which is nearly constant at approximately 12°C. The same well source was used to feed both ponds to control for nutrient quantity and quality. A high-pressure regenerative blower was located within each pond to control for aeration and circulation. Each pond contained 300 mussels of each species from the same age class separated evenly into three suspended upwelling bucket culture systems (SUPSYS) used in standard culture practices in a mixed species assemblage. All mussels of each species used in the chronic temperature experiment were from the same cohort and mussels were randomly collected at each time point. At the start of the exposure (time 0), V. lienosa were of an average length of 25.49±1.75 mm, while V. nebulosa were of an average length of 19.59±1.75 mm. Temperature readings were taken every hour using HOBO^® RH Temperature Loggers (Forest Suppliers, Inc.) submerged to the depth of the buckets.

Mussels were sampled every 2 weeks for 112 days (∼4 months) from July to October 2012. To prevent any variation due to fluctuations from diurnal cycles, all sampling began at 10:00 h and was completed within 90 min. At each sampling time point, mortality was determined for each species per treatment and analyzed as percent cumulative mortality over time. Growth rates were measured as dry tissue and dry shell mass in grams over time for six individuals per species per treatment. These same individuals were then immediately used for determination of overall condition using the following equation: condition index=(dry tissue mass/dry shell mass)×100 (Beninger and Lucas, 1984).

A second subset of six mussels per species per treatment was also sampled every 2 weeks. The foot tissue was dissected for biochemical assays (glycogen and triglycerides) and visceral mass tissue was utilized for gene expression profiling. Foot tissue was immediately flash frozen with dry ice and ethanol after dissection and stored at −80°C. Visceral mass tissue was stored in RNAlater (Qiagen) overnight at 4°C then transferred to −80°C until RNA isolation. Glycogen and triglycerides in the individual foot tissue were measured spectrophotometrically using the Glycogen Assay Kit (Abcam) and Triglyceride Colorimetric Assay Kit (Caymen Chemical Company) following the manufacturer's instructions.

HSP90, HSP70, HSP60 and phosphofructokinase (PFK) were selected for gene expression profiling from the previously constructed Illumina transcriptomic database (Samantha L. Payton, Differential response mechanisms to acute and long-term simulations of global warming for two closely related species of unionid freshwater mussel, Villosa lienosa and Villosa nebulosa, PhD thesis, University of Alabama, 2016). Each transcript was cloned and sequenced for the full-length copy, both to improve primer design for qRT-PCR and to confirm the gene annotation from BLASTX. Based on high sequence homology, a single pair of primers for each gene was designed for the two species (Table 1). Gene expression profiling was performed in the same manner as described for the acute heat shock with the exception that genes were normalized to β-actin after verification that there were no significant treatment effects on β-actin expression. Gene expression profiles were generated from six individual mussel RNA samples per species per treatment instead of pooled replicates.

Statistical analysis of field pond data was performed using SAS Software (version 9.3) and graphs were created with Prism 5 software (GraphPad Software Inc.). To compare the response of each species to the field pond simulations, an analysis of covariance (ANCOVA) with a Tukey's adjustment and Kenward–Roger approximation was performed using PROC MIXED in a repeated measures AR(1) mixed-model on the combined data from all time points and again on the data separated into the heat shock and recovery period (see Fig. 4C for statistical delineation of heat shock and recovery period data). Species, time and treatment were labeled as fixed effects, time was further treated as a repeated measure, and both SUPSYS as well as SUPSYS×treatment were treated as random effects. Variables were transformed when necessary to meet the assumptions of normality and improve homoscedastic model residuals. If variables were transformed in the model, they are represented as back-transformed means and standard errors. For the ANCOVA, significance was based on α=0.10 rather than the usual 0.05 level to allow for all possible significant effects (Burnham and Anderson, 2002; Johnson, 1999). However, it should be noted that the majority of observed significant effects did have P<0.05 (see Tables S2–S6 for full ANCOVA/LSMEANS tables).

An acute temperature challenge was performed to determine whether there were significant differences in the general heat shock response between the Villosa species. Both species were held at either ambient (22°C) or elevated (30°C) temperatures and sampled at 3, 6 and 24 h during the challenge and after a 48 h recovery period (22°C). RT-qPCR was used to determine changes in expression of three predicted heat shock protein genes (HSP90, HSP70-1, HSP70-2). Determination of changes in HSP70 and HSP90 gene or protein expression is routinely used for the assessment of thermal stress based on their important roles in maintaining cellular homeostasis and conveyance of thermal tolerance (Chapple et al., 1998; Sagarin and Somero, 2006; Place et al., 2008; Anestis et al., 2008; Tomanek, 2010). HSP70 proteins, both constitutive and inducible isoforms, are well known for their role in stabilization and assistance in protein folding as part of a normal housekeeping function as well as a protective role during times of cellular stress (reviewed in Daugaard et al., 2007). HSP90 is most notable for its role as a chaperone for mediators of signal transduction pathways, including many involved in cellular stress pathways such as nuclear factor-κB (NF-κB) (Broemer et al., 2004), nuclear steroid receptors (reviewed in Echeverria and Picard, 2010), the aryl hydrocarbon receptor (AHR) (Bell and Poland, 2000) and heat shock factor 1 (HSF-1) (Hu and Mivechi, 2003).

All three heat shock protein genes displayed significant induction in response to the temperature challenge, peaking in expression level at the 6 h time point (Fig. 1). HSP90 expression was significantly upregulated in both species at the 3 and 6 h time points, whereas both HSP70 isoforms were significantly upregulated in both species throughout the 24 h temperature challenge. After 48 h of recovery, HSP90 transcript levels for both species were not significantly different between ambient and warmed treatments. However, V. lienosa exhibited significantly increased expression of both HSP70 isoforms after 48 h of recovery (Fig. 2).

For a more in-depth investigation of the differential expression of the HSP genes between the two species, we compared basal transcript levels at ambient temperature before the start of the temperature challenge. Basal expression of HSP90 was significantly higher in V. nebulosa than in V. lienosa (Fig. 3A). In contrast, V. lienosa had significantly higher transcript levels of HSP70-1, but there was no difference in HSP70-2 transcript levels between the species (Fig. 3A). To assess species-specific differences in the induction of these genes during the acute heat shock, we calculated the ratio of the number of transcripts during heat shock to that under ambient conditions. There was a significant effect of time and species on the induction response for the HSP90 and HSP70-1 genes; however, there was only a significant effect of time for the HSP70-2 gene (see Table S1 for full ANOVA results). There was no significant species-specific difference in induction at the 3 or 24 h time point for any gene tested. However, V. lienosa showed significantly greater induction of the HSP70-1 gene at the 6 h interval (Fig. 3B). Furthermore, V. lienosa had significantly greater induction of HSP90 and HSP70-2 after 48 h of recovery (Fig. 3B).

To determine whether predicted global warming scenarios for the southeastern USA could impact extant populations of freshwater mussels, a chronic warming simulation in field conditions was performed using two culture ponds located at the AABC. One pond simulated standard ambient temperatures and conditions of ponds at the aquaculture facility while a second pond was manipulated in order to create a targeted warmer environment of 3°C above the ambient pond temperature. The warmed pond was successfully maintained at an average temperature 2.5°C above that of the ambient pond for approximately the first two and a half months (July to mid-September) covering the hottest period of the summer (Fig. 4A). By mid-September (after day 70), both pond temperatures gradually decreased with cooling air temperatures and converged between day 98 and day 112. As a result of this convergence of temperature, we were able to analyze the subsequent datasets for either the overall challenge period (day 0–112), or as separate heat shock (day 0–70) and recovery (day 84–112) periods. This proved particularly insightful as the changes in the molecular, metabolic and physiological parameters observed during the overall challenge period were largely driven by the heat shock period.

During the 4 month period of the chronic heat shock experiment, V. lienosa did not experience dramatic mortality in either the ambient or warmed pond (Fig. 4B). However, while V. nebulosa did not experience significant mortality in the ambient pond, there was significant mortality within the warmed pond by day 28 (13% mortality), culminating in 32% mortality by day 70 (Fig. 4C). This mortality correlates to the hottest period within the experimental pond challenge. After 70 days, there was no significant mortality for either species in the ambient and warmed ponds.

To understand how predicted global warming might differentially impact thermally tolerant or sensitive species of freshwater mussel, ANCOVA were performed on data for each growth metric (tissue growth rate, shell growth rate and condition index) for the entire challenge period (day 0–112), and with the challenge divided into a heat shock period (day 0–70) and a recovery period (day 84–112) (see Tables S2–S6 for full ANCOVA/LSMEANS tables). For the combined dataset, the warming treatment did not significantly affect condition or growth rates for either species (Fig. 5). However, growth rates were consistently higher for V. lienosa compared with V. nebulosa within each treatment for the combined dataset and during the heat shock period (day 0–70). During the heat shock period, condition index significantly increased for V. nebulosa in the warmed treatment while there was no significant difference for V. lienosa. During the recovery period (day 84–112), condition index did not differ between the ambient and warmed conditions for either species. Growth rates did not show species-specific trends within treatments during the recovery period, although growth rates for V. nebulosa showed significant declines in the warmed treatment compared with the ambient treatment during the recovery period.

To assess species-specific responses to chronic warming simulation in terms of total glycogen and triglyceride levels, ANCOVA were performed on the data from each tissue biochemistry metric for the entire challenge period (day 0–112), and with the challenge divided into a heat shock (day 0–70) and recovery period (day 84–112) (see Tables S2–S6 for full ANCOVA/LSMEANS tables). For the combined dataset, glycogen levels for V. lienosa were significantly increased in the warmed treatment compared with ambient conditions (Fig. 6A). Although there was no difference in glycogen content between the two species within the ambient treatment, V. lienosa had significantly higher levels within the warmed treatment compared with V. nebulosa. There was a similar pattern during the heat shock period, where V. lienosa possessed significantly higher levels of glycogen within the warmed treatment compared with the ambient treatment, though it was not accompanied by species differences within treatments. The recovery period showed no significant species- or treatment-level differences. For combined data, the warming treatment did not significantly affect triglyceride content in either species, though V. lienosa had consistently higher triglyceride levels in both treatments compared with V. nebulosa (Fig. 6B). During the heat shock period there were higher triglyceride levels for V. lienosa within both treatments compared with V. nebulosa. However, V. nebulosa also showed significantly higher triglyceride levels within the warmed compared with the ambient treatment. The recovery period showed no significant species- or treatment-level differences.

To assess possible species-specific effects on changes in gene expression from the chronic warming simulation, ANCOVA were performed on the data for each molecular biomarker for the entire challenge period (day 0–112), and with the challenge divided into a heat shock (day 0–70) and recovery period (day 84–112) (See Tables S2–S6 for full ANCOVA/LSMEANS tables). In addition to HSP70 and HSP90 expression, we also chose to assess changes in expression of PFK and mitochondrial HSP60 during the chronic challenge. HSP60 plays a prominent role in mitochondrial homeostasis and bioenergetics by assisting in the translocation and proper folding of mitochondrial proteins (reviewed in Deocaris et al., 2006) and participating in the mitochondrial unfolded protein response in times of cellular stress (Yoneda et al., 2004). PFK was chosen for its classic role as a major rate-limiting enzyme in the glycolytic pathway.

For the combined data, the warming treatment caused a significant reduction in the abundance of HSP90 transcripts for V. nebulosa, but did not affect the transcript levels for V. lienosa (Fig. 7A). Also, transcript levels were much higher for V. lienosa in both treatments. During the heat shock period, HSP90 transcripts were significantly reduced for V. nebulosa in the warmed treatment compared with the ambient treatment but were significantly upregulated in V. lienosa. The recovery period showed no significant treatment-level effects; however, V. lienosa had significantly higher HSP90 transcript levels in the ambient treatment. For the combined data, the warming treatment caused a significant decline in HSP70 levels for V. nebulosa and a significant increase in HSP70 levels for V. lienosa compared with the ambient treatment (Fig. 7B). While transcript levels were similar between species in the ambient treatment, V. lienosa had significantly more HSP70 compared with V. nebulosa in the warmed treatment. Although this pattern was similar during the heat shock period, V. nebulosa expressed significantly more HSP70 transcripts compared with V. lienosa in the ambient treatment. While the recovery period yielded no significant treatment differences for either species, V. lienosa had significantly more HSP70 transcripts in both treatments compared with V. nebulosa. For the combined data, the warming treatment affected HSP60 levels of both species compared with ambient treatment levels, resulting in increased expression for V. nebulosa and a decreased expression in V. lienosa in the warmed treatment (Fig. 7C). Although there was no difference in transcript levels between species in the warmed treatment, V. lienosa had significantly higher levels of HSP60 compared with V. nebulosa in the ambient treatment. During the heat shock period, there was still a significant upregulation of HSP60 for V. nebulosa in the warmed treatment compared with the ambient treatment, but there were no significant treatment effects for V. nebulosa. While the recovery period showed no species-specific differences in transcript levels, V. lienosa experienced a downregulation of HSP60 transcript levels within the warmed treatment compared with the ambient treatment. For the combined data, warming had no significant impact on PFK levels for either species, although V. lienosa did have higher transcript levels compared with V. nebulosa in both treatments (Fig. 7D). Although this pattern was similar for the heat shock period, there was no observed treatment or species-specific effect on PFK transcript levels during the recovery period.

Comparative studies related to climate change issues have been largely successful in identifying traits and response patterns useful for informing conservation and management strategies (Klein et al., 2014; Stillman and Somero, 2000; Worthington et al., 2014). Although rather common in marine bivalves, there has been very little comparative research on thermal stress for freshwater mussels. Furthermore, what is available has focused on acute, laboratory-controlled thermal challenges at extreme temperatures (Wang et al., 2012; Luo et al., 2014; Buschini et al., 2003; Pandolfo et al., 2010). Therefore, the main objective of this study was to compare the response of a hypothesized thermally tolerant species (V. lienosa) with a thermally sensitive species (V. nebulosa) of unionid freshwater mussel to both an acute, laboratory-controlled thermal challenge and a field chronic warming simulation at environmentally relevant temperatures. Furthermore, we aimed to illuminate potential differential response patterns between two species that could prove useful for the conservation and management of imperiled freshwater mussel species such as V. nebulosa.

As a first step to understanding how these two species of freshwater mussel might differentially respond to future thermal environments, we subjected them to 24 h acute heat shock at 6°C above ambient laboratory temperatures, coupled with a 48 h recovery period. Because of the short length of this thermal challenge, we focused solely on the molecular response by assessing the expression of classic heat shock protein genes from both species given their well-documented role as a first line of defense against cellular damage (Lindquist and Craig, 1988) and significant role in organismal adaptation to thermal environments (Somero, 2010; Feder and Hofmann, 1999). These factors play an important role in developing acute response models related to gene expression patterns (Wu, 1995; Morimoto, 1998). Specifically, thermally tolerant species may possess higher basal transcript levels and a less pronounced induction profile of heat shock proteins compared with thermally sensitive species (Gleason and Burton, 2015). This could allow thermally tolerant species to withstand a strong thermal insult because of the greater standing pool of heat shock proteins resulting from higher transcript levels (Bedulina et al., 2013; Barshis et al., 2013; Fangue et al., 2006; Luo et al., 2014).

Of the two HSP70 transcripts assessed in this experiment, only HSP70-1 had higher basal expression in V. lienosa (the proposed thermally tolerant species), whereas HSP70-2 levels were equivalent between the two mussel species. Sequencing and complete annotation of both Villosa genomes with additional assessment of HSP70 gene expression profiles would be required to make conclusions about the potential role of differential expression and gene copy number in mediating adaptation to thermal stress in these species. However, some organisms may possess higher copy numbers of key stress protein genes, like HSP70, as an adaptive mechanism for survival in highly variable environments. As a result of high copy number, individual HSP genes could develop new functions and become expressed in response to a stressor or specialized expression under certain conditions (Näsvall et al., 2012). For example, the Pacific oyster (Crassostrea gigas) genome has 88 copies of the HSP70 gene, each with highly variable gene expression profiles in response to cellular stress (Zhang et al., 2012). There is also evidence for species-specific differences in heat shock protein gene copy numbers related to thermal niche in Drosophila species (Evgen'ev et al., 2004).

Basal levels of HSP90 were lower for V. lienosa, deviating from the proposed model expectations of acute response to thermal stress. Most studies comparing HSP90 responses primarily focus on within-species population comparisons (Farcy et al., 2009; Fangue et al., 2006) rather than between-species comparisons. Thus, our understanding of how well HSP90 would follow a traditional heat shock response model in wild populations of thermally sensitive or tolerant species requires further investigation, especially given the complicated nature of the heat shock protein response (as reviewed by Morris et al., 2013). Differential expression of various HSP genes could be responding to trade-offs between maintaining specific HSPs at high basal levels for protection versus the associated energetic costs of maintaining those levels. Such trade-offs would affect the quantity of individual HSP isoforms produced in response to a thermal stress event, and could be due to locally adapted habitat preferences (Tomanek, 2010; Bedulina et al., 2010).

Villosa lienosa has been demonstrated to be tolerant to extended drought conditions (Golladay et al., 2004). Although there are no previous data regarding differential sensitivity to thermal stress between the two species, significant mortality observed in V. nebulosa during the chronic thermal challenge supports the hypothesis that V. nebulosa is more thermally sensitive. This is especially relevant as an increase of only ∼3°C above current water temperatures was enough to cause 33% mortality, suggesting that V. nebulosa may be nearing its upper thermal tolerance limits during peak summer months in its natural habitat. Even though heat shock was more pronounced in the acute heat challenges, mortality differences were not observed. This highlights the importance of using chronic challenges to fully understand how various stressors impact organismal performance.

Our chronic warming simulation was performed at predicted IPCC warming scenarios for the southeastern USA for the year 2100. The resulting data strongly support the conclusion that a predicted thermal increase of 2–3°C is enough to cause differential species-specific responses in physiology, tissue biochemistry and gene expression that could be based on thermal tolerance thresholds. Physiologically, the thermally sensitive species, V. nebulosa, displayed increased condition indices during the heat challenge portion of the simulation. However, though high condition indices indicate rapid increases in tissue mass relative to shell growth, this was not accompanied by significant increases in growth rates or glycogen stores during the heat challenge. The only concurrent observed response was a significant increase in triglyceride stores. It is well documented that upper thermal tolerance limits induce physiological and behavioral modifications that initially lead to higher metabolic rates due to increased oxygen demand (Braby and Somero, 2006; Klepsatel et al., 2013). However, maintaining these high metabolic rates is costly and energetic reserves required for aerobic respiration are eventually depleted (Axenov-Gribanov et al., 2014; Grieshaber et al., 1994), leading to metabolic depression and a shift from aerobic to anaerobic respiration (Anestis et al., 2007, 2008; Sokolova et al., 2012). This could explain the absence of glycogen accumulation, a critical energy reserve for bivalves, and subsequent decrease in growth rates during the recovery period. Conversely, thermally tolerant V. lienosa was completely unresponsive for all physiological parameters, indicative of higher thermal limits required to elicit measurable responses. However, unlike V. nebulosa, there was a significant increase in glycogen reserves during the heat challenge portion of the simulation with no concurrent change in triglyceride concentrations for V. lienosa. This build up of glycogen could represent a decline in metabolism, as glycogen is not being broken down via glycolysis. This glycogen could then be allocated toward maintaining reproductive function over growth or performance (Guderley, 2004; Zani et al., 2012). This is a clear advantage as these two species of unionid freshwater mussel are long-term brooders of developing glochidia from autumn to spring, and long-term environmental stress such as temperature and oxygen limitations could lead to premature evacuation of glochidia (Aldridge and McIvor, 2003).

On a molecular level, there were consistent species-specific patterns in response to the chronic warming simulations. The thermally sensitive species, V. nebulosa, had decreased expression of large molecular weight heat shock proteins HSP90 and HSP70, while V. lienosa experienced an increase in expression for both genes during the chronic heat shock. This is in direct contrast to the acute thermal challenge in which both species displayed significant induction of HSP70 and HSP90. The smaller heat shock protein HSP60 showed the opposite pattern, being upregulated in response to warming for V. nebulosa but showing no response for V. lienosa. As for many of the observations during the chronic thermal challenge, the significant increase in HSP60 expression in V. nebulosa was largely driven by the heat shock period (days 0–70) and there was no significant difference between the ambient and warmed treatments during the recovery treatment. The mitochondrial HSP60 isoform utilized in this study is a biomarker of mitochondrial stress and increased respiratory demands (Hood et al., 2003; Sanyal et al., 1995; Chow et al., 2012), which correlates with the physiological and biochemical results previously discussed.

Previous studies have demonstrated that PFK and other glycolytic or gluconeogenic enzymes, such as pyruvate kinase (PK) and phosphoenolpyruvate carboxykinase (PEPCK), respectively, can be regulated by differential phosphorylation status (Biethinger et al., 1991; Canesi et al., 1999). However, PK and PEPCK have also been shown to be regulated at the transcriptional level in response to stress (Le Moullac et al., 2007). Long-term thermal challenges can lead to metabolic depression, which results in a decrease in aerobic respiration and a greater reliance on anaerobic pathways. Typically, one of the best biomarkers for such an observation is a decrease in the PK/PEPCK ratio (Sokolova and Pörtner, 2001; Múgica et al., 2015). However, PFK is a major rate-limiting enzyme in the glycolytic pathway and may also be regulated at the level of transcription in mussels to meet changing metabolic needs. Villosa lienosa had consistently higher PFK expression levels compared with V. nebulosi, suggesting a higher metabolic demand that is supported by the overall greater growth. Although mollusks may undergo metabolic depression to switch to more anaerobic metabolism, they can also save energy and preserve substrates by lowering their metabolic rate while still remaining in an aerobic metabolic state (Ortmann and Grieshaber, 2003). Although there was no statistically significant change in PFK expression during the heat shock period, the trend toward lower expression may be indicative of a reduced metabolic rate. Additional measurements would be required to determine whether there was any switch from aerobic to anaerobic metabolism.

In general, based on each species' performance during the chronic thermal challenge, the thermally tolerant V. lienosa was able to mount stronger responses without utilizing critical resources necessary to maintain other cellular, metabolic and physiological processes. The thermally sensitive V. nebulosa appeared to have a greater need to utilize these resources in response to thermal stress, although it was ultimately a weaker response. Villosa nebulosa was not able to maintain physiological performance and suffered significant mortality, likely due to a depletion of resources. In contrast, V. lienosa had the resources to maintain a sustained response during the thermal stress and recovery periods. This recovery response is indicative of a thermally tolerant species exhibiting a pre-adaptive or acclimation response (Farcy et al., 2009; Brun et al., 2009; Clegg et al., 1998; Hamdoun et al., 2003), which a thermally sensitive species might not be able to produce. This is an important distinction as acclimation to future thermal warming could help buffer some of the negative impacts for thermally sensitive species (Stillman, 2003). However, to date we have only characterized the thermal stress response of V. lienosa and V. nebulosa from single populations within the state of Alabama. Future studies will investigate responses to thermal stress from multiple populations along longitudinal and latitudinal habitat gradients to address whether the thermal tolerance observed in V. lienosa is more related to greater individual physiological plasticity or greater local adaptation and specialization of genetically structured subpopulations. We will also further investigate the specific mechanisms responsible for any differences in the physiological robustness between the two species that may account for limitations in range expansion by V. nebulosa.

Climate change will continue to impact ecosystems regardless of whether actions are taken to ameliorate some of its causes (IPCC, 2014). Drastic changes to the Earth's terrestrial, marine and freshwater ecosystems have already become evident and are richly reported in the literature. It is crucial to understand how the world's biodiversity is going to respond to future environmental conditions, especially for sentinel and keystone species such as unionid freshwater mussels. Because of long life spans, sedentary lifestyles and obligatory parasitic life stages, migration and adaption might not be viable options for future persistence. Our data reveal that closely related species with overlapping ranges do show species-specific responses to the same environmentally relevant temperature challenges. Further, our data show clear disparities in physiological, biochemical and molecular responses that correlate to unique mechanisms of dealing with thermal stress, which could directly influence reproductive fitness and sustainability of populations in the near future. Finally, traditional and time-consuming metrics, such as condition index, might yield false conclusions as to the health of a wild population if the response mechanism is unknown. In summary, long- term studies using multiple metrics are needed to negate the possibility of potentially ineffective and expensive conservation and management strategies.

We would like to thank undergraduates Greg Ellis and Danielle Porter for help with field work and laboratory assays, as well as Michael Buntin from the Alabama Aquatic Biodiversity Center for help with the construction, implementation and maintenance of the experimental and ambient ponds. We would also like to thank Dr Christina Staudhammer (Department of Biological Sciences, University of Alabama) for aiding in the development of the SAS code for the mixed model ANCOVA.