Feeling the heat: variation in thermal sensitivity within and among populations

Climate change impacts all levels of biological organization, from enzymatic processes to ecological interactions, and the ability of individual species to continue to thrive is often linked to thermal physiological performance (Chown et al., 2010; Deutsch et al., 2015, 2020; Pörtner and Farrell, 2008; Somero, 2012; Sunday et al., 2012). Yet, we lack a thorough understanding of physiological performance variation among individuals and populations that may allow them to withstand increasing temperatures. Variation in thermal physiology is due to both plastic and heritable responses on three time scales: (1) acute – an immediate response without active physiological mitigation; (2) acclimatory – a time-dependent response mitigating the acute response; and (3) adaptive – a population response where there is selection for individuals that are less affected by temperature changes. At an individual level, the response to temperature is described as thermal sensitivity or Q₁₀: the reaction rate fold-change of a physiological trait for every 10°C temperature change (Hochachka and Somero, 2002). This may include differences in acute temperature response, where an individual is acclimated to a single temperature and acutely exposed (acute Q₁₀), or differences in temperature response when individuals are acclimated and assayed at each temperature (acclimated Q₁₀) (Havird et al., 2020). Thus, biochemical reaction rates result in an expected acute Q₁₀ response of 2 (Hochachka and Somero, 2002); however, physiological acclimation or evolutionary adaptation can change this Q₁₀ response, either mitigating the acute effects to maintain homeostasis or intensifying it (Bullock, 1955; Gerken et al., 2015; Klein and Prosser, 1985; Leroi et al., 1994; Sokolova and Pörtner, 2003). Overall, these adaptations, involving both biochemical modifications and physiological processes driving acclimation, render populations more fit for their local thermal environments (Crawford and Powers, 1989; Crawford et al., 2020; Eanes, 1999; Graves and Somero, 1982; Hochachka and Somero, 2002; Powers et al., 1993; Somero, 1995). Thus, in general, ectotherms often evolve mechanisms to maintain similar physiological traits among different thermal environments (Addo-Bediako et al., 2002; Conover and Present, 1990; Conover and Schultz, 1995; Crawford et al., 2020; Dayan et al., 2015; Hochachka and Somero, 2002; Pierce and Crawford, 1997b; Schulte, 2015; Somero, 1978).

Common temperature adaptions to mitigate acute temperature effects suggest strong selection that should reduce interindividual variation within populations. Yet, partitioning the variation within and among populations is challenging, requiring many individuals in several populations with studies that take into account all three response time scales (Havird et al., 2020). Without these data, it is unclear whether individuals have similar thermal sensitivity with little interindividual variation or whether multiple physiological responses drive individual variation in thermal sensitivity. While thermal sensitivity variation in some invertebrates has been observed (Leiva et al., 2018; Nespolo et al., 2003), few studies have examined how it varies among individuals and populations or whether thermal sensitivity is consistent across physiological traits (although see Drown et al., 2021). To examine the variation in thermal sensitivity, we examined three questions: (1) does thermal sensitivity vary among individuals within populations?; (2) is interindividual variation in thermal sensitivity shared across traits?; and (3) does thermal sensitivity vary among populations from different habitat temperatures?

To address these questions, we quantified six physiological traits in Fundulus heteroclitus acclimated and measured at 12 and 28°C. Fundulus heteroclitus is a small estuarine teleost fish, widely distributed along the North American east coast (Burnett et al., 2007). The intertidal salt marshes where these fish reside have large daily temperature, salinity and oxygen fluctuations, contributing to F. heteroclitus’ exceptional physiological plasticity and tolerance (Burnett et al., 2007; Crawford et al., 2020; Smith and Able, 2003). Their wide geographic range and tolerance have made this species a model for examining temperature-driven physiological, biochemical and genetic divergence (Burnett et al., 2007; Crawford et al., 2020). Several studies have previously examined temperature adaptation between extreme northern (e.g. Maine, USA) and southern (e.g. Georgia, USA) F. heteroclitus populations. These population extremes exhibit metabolic and thermal tolerance differences (Fangue et al., 2006, 2009a,b; Healy and Schulte, 2012a, 2019; Healy et al., 2018; Oleksiak et al., 2005; Pierce and Crawford, 1997a; Podrabsky et al., 2000). Yet, there is little information on the interindividual variation in thermal sensitivity within and among populations.

Here, we examined thermal sensitivity within and among F. heteroclitus populations using 350 individuals collected from nine populations from New Jersey to Maine, USA that experience up to a 14.9°C difference in habitat temperature (Fig. 1, Table S1). These populations were chosen because of their non-clinal temperature variation, such that geographically close populations have 2.5–4.4°C habitat temperature differences. To examine whether thermal sensitivity is shared among physiological traits, we measured six traits that are likely evolving in response to local environmental temperature (Burton et al., 2011; Healy and Schulte, 2012a; Oleksiak et al., 2005; Pörtner, 2012; Schulte, 2015): whole-animal metabolic rate (WAM); critical thermal maximum (CT_max); and substrate-specific cardiac metabolism (CM) for glucose, fatty acid (FA), lactate–ketone–ethanol (LKA) and endogenous substrates. Thermal sensitivity was measured as acclimated Q₁₀, such that individuals were acclimated and measured at each temperature. Two of these traits, WAM and CT_max, were measured in each individual at 12 and 28°C. Four other traits, substrate-specific CM for glucose, FA, LKA and endogenous substrates, were terminal determinations and thus measured at either 12 or 28°C. This comprehensive dataset allowed us to investigate the temperature response due to both physiological acclimation and evolved differences among populations. We found that the traits covaried with habitat temperature, and this covariation was dependent on acclimation temperature. We also found very large thermal sensitivity variation with almost 10-fold differences among individuals and significant differences in Q₁₀ related to habitat temperature.

Adult Fundulus heteroclitus (Linnaeus 1766) were collected in September 2018 from nine populations using minnow traps and totaled over 350 individuals from New Jersey (NJ), Massachusetts (MA) and Maine (ME). These sites had significant temperature differences despite their geographic proximity (Fig. 1 Table S1). NJ sites included two reference populations surrounding a site impacted by the thermal effluence from a nuclear power station, reported in a previous publication examining physiological variation between these thermal effluence and reference sites (Drown et al., 2021). Raw data from the temperature and physiological measurements from that study were reanalyzed here with data from six other sites to assess physiological and acclimation differences across a wider temperature range. All fish were tagged with unique visual implant elastomer (VIE) tags for identification throughout the study.

HOBO data loggers collected temperature data throughout August 2018. To ensure HOBO data loggers were accurately collecting temperature data in the species' habitat, loggers were submerged in the shallow marsh. However, as a result of tidal fluctuations, loggers were occasionally exposed at low tide, such that they were recording air temperature, rather than water temperature; thus, only high-tide temperature data were used. Additionally, high tide is typically when most fish are present in the marsh, and move in with the tide (Butner and Brattstrom, 1960). The minimum temperature coinciding with daily high tide time was identified using HOBOware software and then averaged across the collection period. These data were used to determine the mean temperature per location, which was analyzed by one-way ANOVA to determine thermal differences between locations.

All individuals were acclimated to 20°C and 15 ppt salinity for 3 months (12 h:12 h light:dark cycle), a 10°C ‘winter’ period for 6 weeks (8 h:16 h light:dark cycle), then to their experimental acclimation temperature at either 12 or 28°C (16 h:8 h light:dark cycle) for at least 6 weeks prior to physiological determinations. Following initial measurements, fish originally acclimated to 12°C were acclimated to 28°C and vice versa for at least 4 weeks. All fish were fed once daily to satiation (Otohime EP-1 feed) and were housed in recirculating aquaria (density of less than one fish per liter) at 15 ppt salinity. Aquaria were maintained by de-nitrifying biofilter with weekly water exchanges. Fish were fasted 24 h prior to all physiological determinations. Handling and measurement procedures were approved by the Institutional Animal Care and Use Committee guidelines (Animal Use Protocol No: 16-127-adm04).

Following the initial acclimation, WAM was measured for each individual at their acclimation temperature. After a 1 week recovery period, CT_max was measured. After measurement at the initial acclimation temperature, fish were acclimated to the alternative temperature, and WAM and CT_max were measured again at the second temperature. Thus, both WAM and CT_max were determined at both 12 and 28°C for each individual. After a minimum 1 week recovery period post-CT_max, fish were killed and substrate-specific CM was measured at the second acclimation temperature. Because of (limited) mortality and/or technical issues (i.e. sensor malfunction or failure during metabolic data collection), sample sizes varied somewhat between measurements. Full sample size data per trait and temperature can be found in Table 1.

WAM was measured with a custom high-throughput intermittent flow respirometer. This respirometer measures 20 individuals per night, alternating flush and measurements over approximately 12 h as in Drown et al. (2020). Briefly, fish were placed in individual 0.30 l glass chambers. Chambers were then closed off, and oxygen concentration was monitored to measure organismal oxygen consumption rates. Following a 6 or 12 min measurement period, for 28 and 12°C, respectively, chambers were flushed to bring oxygen levels back to 100% saturation. These measurement–flush cycles were repeated continually overnight. The final oxygen consumption rate was determined as the lowest tenth percentile value of a distribution of all measurements throughout the night to estimate metabolic rate (Drown et al., 2020).

CT_max was measured as a proxy for maximum thermal tolerance. Ten fish per measurement were placed in a 10 gallon (∼38 l) glass aquarium with 15 ppt seawater starting at their acclimation temperature (12 or 28°C). A metal heating rod and circulating pump were placed in the aquarium to heat the water at a constant rate of 0.3°C min⁻¹ as in previous studies with F. heteroclitus (Bulger, 1984; Bulger and Tremaine, 1985; Fangue et al., 2006), and an NST thermometer was placed in the tank to accurately monitor temperature. The 0.3°C change is similar to changes that naturally occur in marsh estuaries with an 8°C change in 1 h (Bulger, 1984). An air stone was placed in the water to maintain normoxic conditions. CT_max was recorded as the temperature at which fish lost equilibrium and no longer exhibited an escape response for five continuous seconds.

CM was measured at each acclimation temperature as oxygen consumed by heart ventricles over time as in DeLiberto et al. (2020 preprint). Briefly, fish were killed by cervical dislocation, and ventricles were removed and immediately placed in Ringer's glucose heparin solution. Ventricular oxygen consumption rates were measured in individual chambers inside a temperature-controlled water bath through a fluorometric oxygen sensor spot and fiber optical cable connected to an oxygen meter (PreSens). Metabolism was measured using four substrate conditions: (1) 5 mmol l⁻¹ glucose, (2) 1 mmol l⁻¹ palmitic acid bound to BSA (FA), (3) 5 mmol l⁻¹ lactate, 5 mmol l⁻¹ hydroxybutyrate (ketones) and 0.1% ethanol (LKA), and (4) no metabolic substrate (endogenous) (Oleksiak et al., 2005). Glycolytic enzyme inhibitors (10 mmol l⁻¹ iodoacetate and 20 mmol l⁻¹ 2-deoxyglucose) were added to all but glucose substrates to inhibit any background glycolytic metabolism. Each ventricle was measured in all four substrates, rotated among the four chambers, in the above order for a total of 6 min. Metabolic rate was taken as the slope of oxygen consumption over time for the final 3 min of measurement. Any background flux, measured before and after each run, was subtracted from final measurements.

Statistical analyses were conducted using R version 4.1.0 (http://www.R-project.org/) and RStudio version 1.4.1717 (https://www.rstudio.com/products/rstudio/). An accompanying script detailing the analysis is available from GitHub (https://github.com/ADeLiberto/Fundulus_Physiology). To understand the relationship between habitat temperature and physiological traits, variance due to mass was removed by using the residuals from a log linear (WAM) or linear regression (CT_max and CM) of the trait and mass, at each acclimation temperature. For CM, residuals were calculated both for body mass and heart mass, at each temperature and substrate combination. Additionally, to account for trait variance due to other variables, a forward and backward stepwise Akaike's information criterion (AIC) was used to determine the best-fit model for CT_max and WAM (Hurvich et al., 1998). To test equal variance between traits measured at 12°C versus 28°C, an F-test of variance was used. To represent acclimation response and calculate Q₁₀, mass residuals were transformed back into trait units. Assuming a fish based on the average mass across each dataset, a constant was calculated using the linear regression from mass. This constant was then added to the residuals. Thermal sensitivity (Q₁₀) was calculated using these mass-corrected trait values. Of note, here we use thermal sensitivity specifically for the acclimated Q₁₀, in which individuals were both acclimated and assayed at each temperature. For WAM and CT_max, Q₁₀ was calculated per individual measured at both temperatures. For CM, as an individual was only measured at one temperature, the mean Q₁₀ per population at 12 and 28°C was used to calculate the substrate-specific Q₁₀. Additionally, correlations among traits were examined using a Pearson's partial correlation analysis per temperature using mass residuals (WAM and CT_max) and heart mass residuals (CM) at each temperature.

Metabolic measurements (WAM and CM) were assayed at the acclimation temperature (12 or 28°C); thus, assay temperatures differed by 16°C. CT_max was an acute determination in response to warming temperature (0.3°C min⁻¹), starting at the acclimation temperature. The same individuals were acclimated and measured at both temperatures for WAM and CT_max, allowing us to investigate individual thermal sensitivity (as acclimated Q₁₀) in these traits. As mass significantly affected all traits (Table S2), all trait comparisons include body mass (WAM and CT_max) or heart mass (CM) as a covariate to remove interindividual trait variation due to mass. Body mass and heart mass were significantly correlated at both temperatures with an R² of 0.453 and 0.662 for 12 and 28°C, respectively (Fig. S1A). Interestingly, while body mass was not significantly different between the two acclimation temperatures (Fig. S1C), heart mass was ∼25% greater at 12°C than at 28°C (P<0.001, Fig. S1B).

Based on a 16°C increase in temperature, an expected Q₁₀ of 2.0 should result in a 3.0-fold increase in WAM and CT_max at 28°C relative to 12°C. Among all individuals, WAM was significantly greater, by approximately 2.8-fold, at 28°C than at 12°C (Q₁₀=1.98; Fig. 2A), similar to the expected Q₁₀ of 2.0. Yet, there was high WAM variation among individuals, with greater variance at 28°C than at 12°C (F-test, P<0.001; Fig. 2A, Table 1). At 28°C, CT_max was also significantly greater (1.2-fold) than at 12°C (F-test, P<0.001, Q₁₀=1.11; Fig. 2B), but variance was significantly higher at 12°C than at 28°C (P<0.001; Table 1). For CM, all substrates except for LKA displayed significant temperature responses (Fig. 2C). Glucose and FA metabolic rates were greater at 28°C (P<0.001; Fig. 2C); however, endogenous metabolism was lower at 28°C (P=0.01; Fig. 2C).

Mean minimum high tide temperatures were significantly different between all populations except the two NJ reference populations (NRNJ and SRNJ) surrounding the site heated by nuclear power plant thermal effluence (TENJ; Fig. 1B). Local populations distributed over less than 20 km experienced up to 4.4°C temperature differences, and across all populations there was up to a 14.9°C difference in habitat temperatures.

For WAM, at 12°C, habitat temperature did not influence metabolic rate among these nine populations (P=0.848; Fig. 3A); in contrast, at 28°C, metabolism negatively regressed with habitat temperature, such that colder populations had a higher metabolic rate (P<0.001; Fig. 3A). To examine additional covariates, the relationship between habitat temperature and metabolic rate was also assessed using an Akaike's information criterion (AIC) best-fit model. The input AIC model incorporated acclimation order (12 to 28°C or 28 to 12°C), acclimation temperature, mass and interactions between them (full best-fit model in Fig. S2A, Table S2). Although acclimation order significantly impacted WAM, including these additional covariates did not alter our conclusions with respect to habitat temperature effects on WAM. Interestingly, for 12°C acclimation, all populations had a higher metabolic rate in the 12 to 28°C acclimation order group (Fig. S2A). This difference was not as prominent at 28°C; however, the coldest populations (PEAME/DIME) had a higher Ṁ_O2 when acclimated to 12°C first compared with the group acclimated to 28°C first.

CT_max positively regressed with habitat temperature at both 12 and 28°C, such that CT_max was greater for warmer populations (Fig. 3B). To rule out any possible interactions, an AIC best-fit model was run, incorporating acclimation temperature, acclimation order, sex, mass, habitat temperature and interactions. Comparable to WAM, while significant, the inclusion of these covariates and interactions did not alter our conclusions with respect to the effect of habitat temperature on CT_max (full best-fit model in Table S2, Fig. S2B). As for WAM, differences in acclimation order were more prominent at 12°C. Specifically, CT_max was lower in cooler populations (ME/MA) acclimated to 12°C first, but NJ populations were about the same regardless of acclimation order (Fig. S2B).

Given the difference between acclimation temperatures in terms of heart mass (P<0.001), but not body mass (P=0.87; Fig. S1), we examined CM variation using both heart and body mass residuals. With heart mass residuals, CM was significantly greater in fish from warmer habitats for glucose at both 12 and 28°C, LKA at 28°C, and FA and endogenous metabolism at 12°C (Fig. 4). However, with body mass residuals, these significant relationships were not present (Fig. S3). Finally, to investigate relationships among traits, we examined partial trait correlations within acclimation temperature. Interestingly, CM and WAM were not significantly correlated at either 12 or 28°C, except for endogenous metabolism at 28°C. In fact, among all of the physiological traits, there were few significant partial correlations at either 12 or 28°C, except among the CM substrates (Fig. 5). Additionally, correlations among CM were substrate and temperature dependent, such that FA metabolism was positively correlated with glucose and LKA metabolism at 28°C, but not at 12°C. In. contrast glucose–LKA and LKA–endogenous CM were positively correlated at both temperatures.

Thermal sensitivity or Q₁₀ (change in trait with every 10°C increase) for each individual could be determined for WAM and CT_max because they were measured and acclimated to both 12 and 28°C. For WAM, the high interindividual variation at both 12 and 28°C (Fig. 6A,B) produced a large range of Q₁₀ (from 0.62 to 5.42). In contrast, for CT_max, nearly all the interindividual variation occurred at 12°C, with little variation at 28°C (Fig. 6A,C), resulting in a lower Q₁₀ range for CT_max (1.06–1.20). Interestingly, Q₁₀ values for both traits were significantly related to habitat temperature (P<0.05; Fig. 6D,E), with colder populations having a greater Q₁₀ than warmer populations. Furthermore, WAM and CT_max thermal sensitivity were negatively correlated (P=0.0018), and thus individuals with higher WAM Q₁₀ had lower CT_maxQ₁₀ (Fig. S4). For CM, individuals were measured at either 12 or 28°C, and thus the mean Q₁₀ value per population was calculated. As above, colder populations displayed a higher Q₁₀ for most substrates; however, this was only significant for FA metabolism (Fig. 7). In contrast, colder populations had a significantly lower Q₁₀ for LKA metabolism (Fig. 7).

Temperature is an important factor affecting species distribution and ecological interactions (Deutsch et al., 2015, 2020; Pörtner, 2002, 2010; Somero, 2011; White et al., 2012). Particularly among ectothermic species, there is a consensus that metabolic rates are adaptively important and are driven by natural selection (Anderson and Gillooly, 2018; Clarke and Johnston, 1999; Peck and Conway, 2000). Fundamentally, metabolism has an optimum rate that organisms attempt to achieve by physiological acclimation, evolutionary adaptation or both (Clarke, 2006; DeLong et al., 2018; White et al., 2012). Here, we focused on the concept that strong selection favors an optimum, and thus there will be a reduction in phenotypic variation leading to little interindividual variation.

Acute higher temperatures increase physiological rates up to a maximum or pejus temperature, and then these processes quickly decline (Fig. 8, solid lines) (Pörtner, 2010). Acclimation or adaptation to different temperatures can shift this response and thus reduce temperature sensitivity (Fig. 8, solid blue to solid red line) (Pörtner, 2010; Schulte, 2015), yet rarely do either of these responses completely compensate for temperature differences (Fig. 8, solid blue to dashed red line). In this study, there was a 16°C difference in temperature, and the predicted acute Q₁₀ of 2.0 would produce a 3.0-fold physiological difference.

Overall, our data demonstrate that the six physiological traits have widely different acclimation responses to temperature. On average, WAM showed little temperature acclimation compensation (Q₁₀≈2.0) between the temperatures measured here, similar to previous observations in this species where Q₁₀≈2.0 between similar temperatures (15 and 30°C; Healy and Schulte, 2012a). In contrast, CT_max had much lower thermal sensitivity and greater acclimation compensation (Q₁₀≈1.1) between acclimation/measurement temperatures of 12 and 28°C, as previously observed across a range of temperatures in this species, where average Q₁₀ in a New Hampshire (northern) population ranged from ∼1.0 to 1.25 (Fangue et al., 2006). We expected a similar scaling of the temperature response for all metabolic rates (WAM and CM), yet this was not the case. Despite high thermal sensitivity in WAM, CM showed surprisingly small or no thermal sensitivity, with Q₁₀ ranging from 0.90 to 1.18, which also contrasts with maximum heart rate measurements in F. heteroclitus that differ quite widely when acclimated to 15°C versus 30°C (Safi et al., 2019).

The small or insignificant difference across 16°C in CM for FA, glucose and LKA suggests a total compensatory acclimation (Q₁₀≈1.0) (Fig. 8). Endogenous CM was significantly higher at 12°C than at 28°C (Fig. 2) and thus the Q₁₀ was less than one. Cardiac metabolic thermal sensitivity may reflect larger endogenous substrate stores at 12°C, which allows for a greater endogenous metabolism at lower temperatures, similar to what was previously observed across only NJ individuals (Drown et al., 2021). The observed larger heart mass at 12°C, resulting in greater endogenous metabolic stores, would support this supposition. An alternative explanation is that for CM, these individuals lie on either side of a thermal performance curve (Schulte et al., 2011); however, this is unlikely given 28°C is within the standard thermal range for this species, and no F. heteroclitus performance curves have a pejus temperature at or below 28°C (Baris et al., 2016b; Chung et al., 2017; Fangue et al., 2008; Healy and Schulte, 2012a; Johnson and Bennett, 1995).

Conceptually, reduced temperature sensitivity with acclimation allows individuals to maintain a constant performance across a range of temperatures, providing an advantage compared with an acute temperature response. Acclimation may be accomplished by increasing enzyme concentrations at lower temperatures through mRNA regulation. For example, lactate dehydrogenase-B (LDH-B) and other proteins increase in concentration at lower temperatures through mRNA regulation (Segal and Crawford, 1994). Clearly, with a Q₁₀ of ∼1.0, CM is insensitive to temperature change, and thus physiological mechanisms eliminate the acute effect of temperature on biochemical reactions rates. These mechanisms are missing for WAM, where the Q₁₀ is ∼2.0, suggesting that the need to integrate across physiological processes to define overall metabolic rate limits the ability to compensate for temperature changes. However, the interindividual variation we observed in WAM Q₁₀ (as discussed below) suggests that, instead, there are biochemical and physiological mechanisms that readily shift the WAM thermal response, but these responses are highly variable, resulting in a mean Q₁₀ of 2.0.

The small but statistically significant CT_max difference between acclimation temperatures reflects the relatively small variation among all individuals (Fig. 2B, Table 1). The thermal sensitivity in CT_max may depend on the specific acclimation temperature. The small difference among individuals acclimated to 12°C versus 28°C (16°C difference) could reflect a CT_max maximum limit or hard ceiling where both the variance and the mean plateau at an upper temperature (Morgan et al., 2020). Indeed, previous work in F. heteroclitus measuring CT_max at multiple acclimation temperatures (Fangue et al., 2006) showed a plateau in CT_max ∼42°C, which is consistent with the average CT_max observed here (41.8°C). Acclimation to two closer temperatures (<16°C difference) may result in a similar difference in CT_max (i.e. ∼6°C), but in turn a lower Q₁₀.

The benefit or direct fitness effect for lower CT_max at 12°C acclimation is not clear because of the highly variable marsh environment temperatures. For F. heteroclitus, marsh temperatures are likely to exceed 35°C (average CT_max at 12°C; Fig. 2B) and cold-acclimated individuals with a lower CT_max may have reduced fitness when subjected to these temperatures. While the thermal tolerance of marine species is an important ecological parameter for determining species range (Sunday et al., 2012), we suggest this reduction in CT_max with acclimation reflects acclimatory effects on other physiological processes. However, the lack of correlation among the traits measured (Fig. 5) suggests these effects may be more complex or tied to another trait not measured here.

In F. heteroclitus, metabolism and the underlying biochemical processes have evolved to counteract temperature effects (reviewed in Crawford et al., 2020). We measured multiple traits in the same individuals among populations inhabiting a mosaic of temperatures ranging from 15 to 32°C. We found that habitat temperature explained a significant amount of variation in several of these traits (Figs 3 and 4). For WAM, we observed higher metabolic rates in the colder populations when acclimated/measured at 28°C, yet no difference among populations from different habitat temperatures when acclimated/measured at 12°C, even when accounting for additional covariates (e.g. acclimation order; Fig. S2, Table S2). Previous studies examining metabolic rate in F. heteroclitus have been inconclusive. In New Hampshire (NH) and Georgia (GA) populations, Ṁ_O2 was significantly greater in the colder populations at 5, 15 and 25°C (Fangue et al., 2009a). However, in acclimated routine Ṁ_O2 measurements, differences between MA and GA populations were only observed at 5°C, across temperatures from 5 to 33°C (Healy and Schulte, 2012a). Furthermore, Ṁ_O2 in NH, NJ and GA populations measured at 15°C showed significantly lower metabolic rate in GA, but no significant difference between NH and NJ (Brennan et al., 2018; Healy et al., 2019). In our study, the significant relationship between WAM and habitat temperature at 28°C (Fig. 3A) was mostly driven by the three most southern (NJ) populations. Overall, it appears that the impact of habitat temperature on metabolic rate is acclimation/measurement temperature specific. Metabolism may also be lower only in the NJ populations because of the historical divide and admixture zone associated with genetic isolation, including different mitochondrial haplotypes, similar to what has been observed with hypoxia tolerance (Brennan et al., 2018; Crawford et al., 2020; Healy et al., 2018).

In contrast, across all nine populations, CT_max at both temperatures was significantly greater in warmer populations (Fig. 3B). This increase in CT_max in warmer habitats was observed even after long-term acclimation to two temperatures and thus is not due to reversible physiological acclimation. Therefore, it is most likely heritable, although we cannot rule out irreversible developmental or transgenerational effects (Cavieres et al., 2019). Thus, the variation in CT_max across habitats likely represents local adaptation. This supports previous work in F. heteroclitus where significant CT_max differences were observed both at the extremes of its range and between closely related populations (Dayan et al., 2015; Fangue et al., 2006; Healy et al., 2018). Yet, it is unclear why lower CT_max would be favored even at lower habitat temperatures. This may reflect a more derived state. Conversely, CT_max may have underlying physiological pathways involved in other processes where higher rates would be unfavorable.

CM (Fig. 4), unlike WAM, was lower at colder habitat temperatures depending on the substrate and acclimation temperature. These results are similar to heart mitochondrial respiration, where colder populations of Fundulus species have lower metabolism (Baris et al., 2016b, 2017). Based on the data presented here, this habitat temperature response is not due to the inability to modify CM as exemplified by the similar cardiac metabolic rate when acclimated and assayed at 12 and 28°C. That is, after acclimation, CM Q₁₀ values are close to 1.0 and therefore show similar rates with a 16°C assay temperature difference; yet, there was a significant decrease with ∼12°C habitat temperature difference. Overall, these six traits show significant divergence among populations related to habitat temperature. These traits have been shown to be heritable (Crawford et al., 2020; Pörtner, 2012; Schulte, 2015), and it is likely that local adaptation is driving phenotypic shifts based upon temperature selection.

Additionally, while individuals were considered to be fully acclimated after both acclimation phases, we did observe an acclimation order effect on both WAM and CT_max, which was also previously observed in only NJ individuals (Drown et al., 2021). While this effect was significant, technical variation in acclimation order was corrected for using an AIC best-fit model for both CT_max and WAM, and habitat temperature remained a significant factor explaining trait variation for all but WAM at 12°C, which was also not previously significant (Fig. 3; Fig. S2).

While local adaptation of physiological processes is common, especially in F. heteroclitus (Crawford et al., 2020), and the within-population variation in these specific physiological traits that we observed has been found in other independent studies (Healy and Schulte, 2012a; Healy et al., 2018; Oleksiak et al., 2005), our large dataset uniquely allowed us to examine variation in temperature response and thermal sensitivity within and among populations for a suite of metabolic and thermal tolerance traits.

To understand the variation in metabolic rates among individuals, we compared individual Q₁₀ values for WAM from the individual data. Surprisingly, Q₁₀ ranged from 0.62 to 5.42. Thus, some individuals had the same or lower WAM rates at 12 and 28°C, while others had a Q₁₀ that exceeded the expected acute temperature effect (Fig. 6). This variation in Q₁₀ is unexpected for two reasons. First, biological processes are expected to evolve to maintain homeostasis and, pertinent to this study, mitigate responses to environmental temperature variation. Second, thermal sensitivity should be adaptive and natural selection for an optimum would reduce individual variation within a population. That is, selection for a specific trait should favor an allele or combination of alleles, thus reducing the frequency of the alternative allele, which would reduce heterozygosity or genetic variation and in turn the heritable trait variation. Yet, we found that some individuals had almost total compensation (nearly equal metabolic rates at 12 and 28°C, Q₁₀≈1.0), while others had an almost 15-fold increase, yielding a Q₁₀ of 5.4, greatly exceeding the expected Q₁₀ of ∼2.0 in metabolism (Fig. 6A). We suggest that this variation among individuals represents different strategies for coping with environmental temperature variation. Some individuals may overcome the physical effect of temperature, while others exploit these higher temperatures.

Thermal sensitivity was significantly correlated with habitat temperature for many of the traits we measured, suggesting that it is both biologically relevant and adaptively important (Fig. 6D,E). The Q₁₀ for both WAM and CT_max across habitat temperatures indicates lower sensitivity among populations from warmer habitats. For WAM, this pattern appears to be driven by the NJ population. While the ME and MA populations appear to have a similar Q₁₀ values, those for the NJ population were much lower. As discussed above, the NJ populations are south of a historical evolutionary break, and this north–south historical isolation may be driving this pattern. In contrast, for CT_max, the pattern is more linear, supporting the idea that local habitat temperature is driving the Q₁₀ response. CM, with the exception of that for LKA (Fig. 7), was lower at lower habitat temperatures depending on the substrate and acclimation temperature, unlike WAM, but similar to heart mitochondrial respiration, where colder populations of Fundulus species have lower metabolism (Baris et al., 2016a, 2017).

Overall, WAM and CT_max measurements at two acclimation and assay temperatures among the same individuals indicate large interindividual variation that results in a wide range of Q₁₀ values. While only a small percentage of this variation can be explained by habitat temperature (Fig. 6D,E), patterns among populations suggest that the Q₁₀ variation is biologically relevant. Yet, most of the Q₁₀ variation is among individuals within a population. We suggest that individual variation in thermal sensitivity results in different physiological strategies in response to environmental temperature variation, which may promote the maintenance of standing genetic variation in a population (Burton et al., 2011; Careau et al., 2014; Norin et al., 2016). Furthermore, our findings regarding lower thermal sensitivity at higher temperatures support previous data suggesting individuals at lower latitudes have lower thermal sensitivity than higher latitude individuals (Seebacher et al., 2015).

Our data specifically represent thermal sensitivity between acclimation and measurement temperatures of 12 and 28°C, which leaves much of the F. heteroclitus thermal range unexplored in the current study. Thus, it is likely that variation in Q₁₀ may depend on the temperature range being examined. For example, in the study by Healy and Schulte (2012a) of metabolic rate across a range of temperatures, thermal sensitivity was much more dependent on the temperatures being examined. Yet, Fangue et al. (2009b) found the thermal sensitivity of CT_max remained somewhat consistent (Q₁₀≈1.1) across a 30°C temperature range. Lastly, repeated measurements in an individual for Q₁₀ would provide insight as to whether this trait is repeatable, although given past repeatability of metabolic rate and CT_max, it is likely to be consistent (Drown et al., 2020; Healy and Schulte, 2012b; Morgan et al., 2018; Nespolo and Franco, 2007).

Global climate change will require species, populations and individuals to adjust to rapidly changing environments. If the physiological responses we measured are heritable, the data presented here on the individual variation in physiological traits and their Q₁₀ values support a large standing genetic variation, particularly if individuals have different physiological strategies to cope with change. This breadth of standing genetic variation would enhance rapid evolution in physiological performance to compensate for global climate change (Matuszewski et al., 2015; Scheffers et al., 2016), providing some hope for species survival.

The authors would like to thank Liam Dorsey, Agatha Freedberg and Rebecca VanArnam for assistance with data collection and fish husbandry, and John Proefrock for comments on the manuscript.