Thermodynamic effects drive countergradient responses in the thermal performance of Littorina saxatilis across latitude

We compared cardiac thermal performance curves (TPCs)


H I G H L I G H T S G R A P H I C A L A B S T R A C T
• We compared cardiac thermal performance curves (TPCs) across latitude. • High-latitude populations possessed greater upper performance and tolerance limits. • Countergradient responses driven by elevated metabolism in cold-adapted gastropods. • Latitudinal differences were less pronounced following warm acclimation. • Thermal selection at low temperatures can have knock-on consequences for upper limits.

A B S T R A C T A R T I C L E I N F O Editor: Rafael Mateo Soria
Keywords: Countergradient Latitudinal gradients Thermal performance curve Thermodynamic effects Metabolic cold adaptation Local adaptation Thermal performance curves (TPCs) provide a powerful framework to assess the evolution of thermal sensitivity in populations exposed to divergent selection regimes across latitude. However, there is a lack of consensus regarding the extent to which physiological adjustments that compensate for latitudinal temperature variation (metabolic cold adaptation; MCA) may alter the shape of TPCs, including potential repercussion on upper thermal limits. To address this, we compared TPCs for cardiac activity in latitudinally-separated populations of the intertidal periwinkle Littorina saxatilis. We applied a non-linear TPC modelling approach to explore how different metrics governing the shape of TPCs varied systematically in response to local adaptation and thermal acclimation. Both critical upper limits, and the temperatures at which cardiac performance was maximised, were higher in the northernmost (cold-adapted) population and displayed a countergradient latitudinal trend which was most pronounced following acclimation to low temperatures. We interpret this response as a knock-on consequence of increased standard metabolic rate in high latitude populations, indicating that physiological compensation associated with MCA may indirectly influence variation in upper thermal limits across latitude. Our study highlights the danger of assuming that variation in any one aspect of the TPC is adaptive without appropriate mechanistic and ecological context. Science of the Total Environment 863 (2023) 160877

Introduction
Latitudinal variation in temperature has a pervasive impact upon the performance and distribution of organisms (Gaston et al., 2009;Hochachka and Somero, 2002;Spicer and Gaston, 1999). Comparisons between populations and species from different latitudes represent an important means to study the adaptive physiological mechanisms enabling organisms to cope with different thermal regimes (Angilletta, 2009;Leung et al., 2021;Somero, 2005), while also providing valuable insights to assess the relative susceptibilities of different organisms and ecosystems to future climate change (Buckley and Huey, 2016;Deutsch et al., 2008;Han et al., 2019;Morley et al., 2022). Intertidal ectotherms represent a widely utilized model for the study of thermal adaptation across latitude, given that they are subject to wide temperature variation through the tidal cycle and possess latitudinal distributions which are believed to be highly constrained by their thermal tolerance limits (Dong et al., 2021;Han et al., 2019;Somero, 2010). Numerous studies have demonstrated that intertidal species from lower latitudes display more elevated upper thermal limits than those from higher latitudes, matching broad-scale differences in habitat temperatures (Compton et al., 2007;Dong et al., 2021;Stillman and Somero, 2000). Intraspecific comparisons, however, have presented a mixed picture, with some studies identifying intraspecific variation in upper limits across latitude (Gaston and Spicer, 1998;Madeira et al., 2012), but others failing to find evidence of consistent patterns (Dong et al., 2015;Gaitán-Espitia et al., 2014;Jupe et al., 2020;Kuo and Sanford, 2009;Lee and Boulding, 2010;Logan et al., 2012). This picture is further complicated by the existence of mosaic patterns of temperature variation between intertidal sites (Helmuth et al., 2006), driven by factors including tidal cycles, local upwellings and topography. These may obscure or even override the effects of broad-scale latitudinal thermal adaptation at local and intermediate spatial scales (Han et al., 2019;Helmuth et al., 2006;Kuo and Sanford, 2009). Despite the lack of consensus regarding the relationship between latitude and intraspecific local thermal adaptation, meta-analyses comparing the upper thermal limits of different species across latitude frequently assume a single invariant response across a species' latitudinal range (Bennett et al., 2021;Deutsch et al., 2008;Sunday et al., 2014), thereby making the implicit assumption that thermal limits within species are highly conserved. Such assumptions risk underestimating the potential role of intraspecific local adaptation in determining how populations may respond differently to thermal stress across latitude, particularly in widely-distributed species (Bennett et al., 2019;Spicer and Gaston, 1999). This, in turn, may have substantial repercussions for our understanding of the potential ecological impacts of future climate change, including both potential impacts on local abundance and species range limits (Bennett et al., 2021;Buckley and Huey, 2016;Han et al., 2019;Morley et al., 2022). Thus, in order to better understand globalscale patterns of thermal adaptation and predict future changes, it is crucial to improve our understanding of the extent and nature of intraspecific variation in upper thermal limits across latitude, as well as its mechanistic underpinnings.
In addition to upper thermal limits, local adaptation of metabolic rate plays a crucial role in responding to different temperature regimes across latitude (Hochachka and Somero, 2002). Organisms from high latitude populations often compensate for the reduction in biochemical reaction rates experienced in colder climates by elevating their metabolism via a process known as 'metabolic cold adaptation' (MCA) (Addo-Bediako et al., 2002;Bozinovic et al., 2011;White et al., 2012). This process may be underpinned by increases in the concentration and/or catalytic activity of enzymes, ultimately increasing metabolic rates relative to warmadapted populations kept at equivalent temperatures (Hochachka and Somero, 2002). MCA is closely linked to countergradient variation, a phenomenon whereby life history traits such as growth, development and feeding rate and body size are elevated in high latitude compared to low latitude populations leading to either a reduction in, or reversal of, patterns of phenotypic variation in these traits across latitude compared to what would be expected based upon the biological effects of temperature (Conover and Schultz, 1995;Conover et al., 2009;Levins, 1968;Levinton, 1983). Compensatory adjustments associated with MCA to cope with low temperatures may also have knock-on consequences for performance across a wide range of temperatures, and so have the potential to lead to shifts in the shape of thermal performance curves (TPCs), both due to trade-offs in resource allocation (Huey and Kingsolver, 1989;Levins, 1968) and thermodynamic effects on metabolism and enzyme stability (Gillooly et al., 2001;Hochachka and Somero, 2002). Crucially, because the temperature at which performance is maximal (T peak ; widely referred to in the literature as the thermal optimum, T opt ) and the critical thermal maximum (CT max ) are integrated components of the TPC (Angilletta et al., 2010), any shifts in the shape of TPCs arising due to MCA could also have important repercussions for latitudinal patterns in upper thermal limits. However, predicting the extent and nature of these repercussions is conceptually challenging, as indirect effects on upper limits could arise via several alternative mechanisms (Angilletta, 2009;Gardiner et al., 2010;Villeneuve et al., 2021). Elevated metabolic rates in high latitude populations could, for instance, be achieved via a leftwards shift in the TPC, leading to reduced performance at high temperatures and lower T peak and CT max (cogradient model) (Angilletta, 2009). Alternatively, elevated metabolic rates in highlatitude populations might be sustained across the entirety of the TPC, potentially leading to greater maximum performance and T peak relative to low-latitude populations (countergradient model) (Angilletta, 2009;Gardiner et al., 2010;Villeneuve et al., 2021). As even small changes in single parameters of the TPC could have important knock-on consequences on upper limits (Rezende and Bozinovic, 2019), it is vital to adopt approaches which allow for all parameters governing the shape of TPCs to be incorporated into an integrated framework.
To advance our mechanistic understanding of how local intraspecific thermal adaptation across latitude may influence the shape of TPCs, and in turn how this may affect upper thermal limits, we assessed systematic variation in thermal performance across four latitudinally distinct populations of the intertidal periwinkle Littorina saxatilis. To distinguish between evolutionary and plastic (acclimation) responses across latitude, we acclimated snails to three different temperatures prior to experiments. We measured cardiac activity in response to an acute thermal ramp, both as a proxy for the thermal performance of metabolism (Frederich and Pörtner, 2000), and as an established method to assess upper thermal limits in intertidal gastropods (Dong et al., 2021). We combined this with a recently developed, non-linear regression approach to fit TPCs (Bozinovic et al., 2020a;Rezende and Bozinovic, 2019), providing us with an integrated framework to assess how multiple performance characteristics vary in tandem to affect the shape of TPCs. L. saxatilis possesses a wide latitudinal distribution (Reid, 1996), direct development coupled with poor dispersal, and high genetic differentiation between populations (Westram et al., 2016). We predicted that the species would display a substantial degree of intraspecific variation in both the shape of TPCs and upper thermal limits, providing us with a suitable model to explore the interrelationship between these traits and thus help elucidate broad-scale patterns of thermal adaptation across latitude.

Collection sites and animal husbandry
We obtained adult populations of L. saxatilis from four intertidal sites (referred to as Iceland, Scotland, England and Spain) across a latitudinal gradient spanning 23°along the northeast Atlantic coastline (Table S1). Sites were chosen with a high degree of latitudinal separation to reduce the effect of "mosaic" patterns of temperature variation, which are known to confound latitudinal patterns of thermal adaptation over intermediate spatial scales in the intertidal. To confirm the existence of a latitudinal trend in temperature, historical air and sea surface temperature data for the Northeastern Atlantic were obtained from the WorldClim and Bio-Oracle online databases, and annual mean air and sea surface temperatures identified for the four locations examined in this study. A clear latitudinal trend can be observed (see supplementary Fig. S1, Table S1, and accompanying captions). Collection dates at different sites were staggered from August to December 2019 to reduce variation in seawater temperatures at time of collection (~13-15°C in the three southern populations). As seawater in Iceland does not reach this temperature, snails from this population were instead collected during the maximum seawater temperatures experienced annually (~8°C). Because microhabitat variation and shore height are known to affect thermal sensitivity in L. saxatilis (Dwane et al., 2021;Sokolova et al., 2000), all collections were made from upper intertidal populations belonging to the "moderate" ecotype sensu Reid (1996), from comparable microhabitats (rocky outcrops with minimal algal and boulder cover). Snails were placed in plastic containers with dampened cardboard, and transported to the University of Plymouth, either by air freight (<72 h) for the Icelandic and Spanish populations, or by car (<24 h) for the English and Scottish populations. Upon arrival, snails were distributed randomly between 5 L aquaria containing aerated sea water (salinity = 35 ± 1, stocking density~10 snails L −1 ). All populations were initially held at 15°C for >1 week to standardize recent thermal history, except for the Icelandic population (10°C). Populations were then acclimated to 10, 15 and 20°C by transferring the aquaria to controlled temperature rooms at the respective temperatures (n ≥ 3 aquaria per treatment) and were held at these temperatures for a minimum of 5 weeks, which is sufficient to limit differences in metabolism arising from acclimatization in the field (Hawkins, 1995;Newell, 1979). The three acclimation temperatures were chosen to encompass a range of habitat temperatures experienced by the species across its latitudinal distribution. Specifically, 10°C approximates the average maximum air and water temperatures experienced by the Icelandic population in summer months (12.20 and 11.08°C, respectively), while being close to the annual mean air and water temperatures experienced by the English (10.4 and 12.27°C) and Scottish (8.4 and 10.95°C ) populations;. 15°C approximates mean air and water temperatures experienced by the Spanish population annually (14.3 and 15.5°C, respectively), while 20°C approximates average maximum summer air temperatures for the three southern populations (23.0, 19.6 and 18.1°C for the Spanish, English and Scottish populations respectively) and average maximum water temperatures for the Spanish population (18.51°C). Water changes were conducted weekly, and snails were fed ad libitum on Ulva lactuca. Snails were able to climb into and out of the water within the aquaria, and thus experienced thermal acclimation to the same temperatures under both submerged and aerial conditions. Aquaria lids were sealed with cling film to prevent escape and maintain humidity. Mortality during the acclimation period was <10 %, except in 20°C-acclimated snails from England (~20 %), 15°C-acclimated Scottish snails (~20 %) and 20°C-acclimated Scottish snails (~45 %).

Ramping protocol
Cardiac performance (as a proxy of metabolic rate) in response to acute thermal ramping was measured in air using non-invasive photoplethysmography (Burnett et al., 2013). Plethysmographic signals were detected using an infrared sensor (CNY70, Vishay semiconductors, Germany) positioned over the heart of the snail and held in place using plastic mesh. This prevented the snails from freely moving during the trial, while not restricting their ability to respire or emerge from their shells. The infrared sensor was connected to an amplifier unit (AMP03, Newshift Lda., Portugal) and data acquisition system (USB-6001, National Instruments, USA) connected to a laptop running NI-DAQExpress (National Instruments, USA). Amplified signals were recorded at a sample frequency of 40 Hz.
After being fitted with sensors, snails (n = 4-5 per trial) were transferred to individual dry plastic bags placed inside a water bath (Grant R5, Grant Instruments, UK) preheated to the required acclimation temperature (10, 15, or 20°C). The bags were weighted with marbles to submerge them within the water bath, but were left unsealed above the waterline to facilitate gas exchange. Preliminary trials confirmed snails experienced normoxic conditions (> 95 % air saturation) throughout the ramp. Following addition of the snails, the setup was left for 15-20 min to allow bag temperatures to equilibrate with water bath temperatures. Temperature during the trials was monitored using a HH806U temperature logger (Omega, USA) with a k-type thermocouple inserted inside an empty snail shell filled with blu-tak (Bostic, UK) and placed inside an additional envelope bag. This was done to mimic body temperature of the snails during the trials by taking into account thermal inertia provided by the snail shells, as well as the air space surrounding the snails within the bags. Following equilibration, temperature was increased at a rate of 6°C h −1 until after all snails reached cardiac flatline (~47-50°C), at which point the experiment was terminated. Shells were then measured (columella height and shell width across the widest part of the shell, perpendicular to the columella). To obtain wet weight, shells were cracked open to dissect the soft tissue, which was rinsed, blotted dry and weighed (± 0.01 mg) using a Cubis Semi-Micro Balance (Sartorius, Germany). Size and weight information for each population is given is supplementary table S1. Three to four trials were run per population per acclimation temperature. Faint or inconsistent heartbeat traces were discarded from further analysis (39 out of a total of 181 individuals used in thermal ramps). This occurred either due to amplifier battery failure, poor positioning of the sensors, or dislodgment of the sensors during the trials.

Data analysis
Raw signals were converted into heartrate data using the MATLABbased Artiifact package (Kaufmann et al., 2011). Heartbeats were counted for 1 min periods at 1°C (± 0.1°C) intervals through the temperature ramp. This was done by detecting voltage peaks in the raw signals using global peak thresholds, followed by manual screening to correct miscounted sections of the heartbeat signal (e.g. due to the presence of doubly-counted peaks, Burnett et al., 2013). Following this, heartrate in BPM was calculated.
We fitted TPCs to heartrate data separately for each individual using the "fit.thermal.curve" function (Bozinovic et al., 2020a;Rezende and Bozinovic, 2019) in R version 4.0.5. This non-linear function, available at Dryad Data Repository (Bozinovic et al., 2020b), fits curves to thermal performance data based on four descriptive parameters: a constant that shifts the curve up and down everything else being equal (C), the fold-increase in performance expected for a 10°C rise in temperature (Q 10 ), the threshold temperature in which denaturation effects are observed (T th ) and the thermal sensitivity of denaturation to temperatures above this value (d) (Rezende and Bozinovic, 2019). Based on these parameters, the descriptors T peak , P max (heartrate at T peak ), and T breadth (the thermal breadth in which performance is a 50 % of its maximum) were obtained from function outputs as described by Rezende and Bozinovic (2019) and Bozinovic et al. (2020a). CT max , was obtained by extrapolating the temperature at which heart rate reached zero from the TPC fits, while P initial was calculated by extrapolating heart rate values at the respective acclimation temperature of each individual (Fig. 1A). In cases where initial convergence could not be achieved for some individuals, we instead employed the "boot.thermal. curve" function (Bozinovic et al., 2020a), also available at Dryad Data Repository (Bozinovic et al., 2020b). This method accounts for slight deviations from the expected TPC shape present in the raw data by iteratively fitting the nonlinear model to bootstrapped data values from a normal distribution ± s.d. around the mean raw data at each temperature, with standard deviation values of 1-5 bpm used depending on the individual in question. The bootstrap analysis was run iteratively 20 times, with median Q 10 , C, T th and decay values used to generate TPCs for these individuals. TPCs were assessed visually for goodness-of-fit to raw data, resulting in the exclusion of four individuals due to unreasonable model fits. In total, successful convergence of the non-linear regression model to heartrate data was achieved across 128 individuals, of which 50 required the use of bootstrapping. Fits of the non-linear model to the raw data were high across both non-bootstrapped and bootstrapped individuals (r 2 = 0.88 ± 0.094, mean ± sd). Plots of the individual non-linear model fits are presented in supplementary fig. S2.
To investigate how widely utilized metrics of thermal tolerance and performance differed in response to latitude and acclimation temperature, as well as to explore the possibility of an interactive effect, General Linear Model (GLM) analyses were performed separately for six standard metrics of thermal performance: CT max , T peak , P max , P initial , Q 10 , and T breadth . These were performed using the base R function "lm" and the "Anova" function in the package car, with type III sum of squares. In all models aside from for T breadth , latitude and acclimation temperature were included as continuous and categorical predictor variables respectively. However, preliminary analysis identified that variation in T breadth followed a markedly non-linear pattern across latitude; therefore, latitude was instead treated as a categorical factor for T breadth . Normality of residuals and homoscedasticity of variance was verified using Quantile -Quantile and scale-location plots using the "plot" function in base R. To meet GLM assumptions, Q 10 was logtransformed prior to analysis. To control for the effects of size differences between populations, wet weight was initially included as a covariate in all five models. However, no significant effect of size was observed for any parameter except P max , thus it was dropped as a covariate prior to final analysis in the other five models. Post-hoc tests were conducted using either the "emmeans" function to compare across acclimation temperatures in the absence of a significant interaction term or using the "emtrends" function to compare slope values in the presence of a significant interaction, both from the R package "emmeans". For T breadth , "emmeans" was used to compare across latitude within each acclimation treatment. In addition, to explore the relationship between the performance traits CT max , T peak and P max across the dataset as a whole, we performed Pearson's correlation tests using the "corr.test" function in base R. To determine how different TPC parameters and descriptors generated by the non-linear model interact to contribute to inter-individual variation in the shape of TPCs within the dataset, we conducted a Principal Component Analysis (PCA) using the R function "prcomp". This incorporated the four parameters Q 10 , C, T th and log-transformed d (to improve normality), as well as the descriptors T peak , CT max , T breadth and P max . It should be noted that while we obtained raw values for CT max based on cardiac flatline during the thermal ramps, we chose to use extrapolated CT max values in both the linear model analysis and PCA for consistency with other values derived from the model fits. Fitted values for CT max were extremely similar to raw values across the dataset (T 126 = 38.93, p < 0.001, Pearson's correlation coefficient = 0.96) and substitution of CT max with the raw values did not qualitatively affect the statistical output of the linear model.

Thermal performance traits
In contrast to the latitudinal trend in habitat temperature (Fig. S1), critical thermal limits (CT max ) displayed a significant positive trend with increasing latitude (F 1,122 = 36.953, p < 0.001; Fig. 2A; see also supplementary Table S2 for full ANOVA outputs), being 1.7°C higher in the northernmost (Icelandic) population than in the southernmost (Spanish) population when averaged across all three acclimation temperatures. There was neither an effect of acclimation temperature (F 2,122 = 0.105, p = 0.9), nor an interactive effect of acclimation temperature x latitude (F 2,122 = 0.027, p = 0.974). There was also a significant positive effect of latitude on T peak (F 1,122 = 5.923, p = 0.016; Fig. 2B) with no significant interaction term (F 2,122 = 1.876, p = 0.157), A similar latitudinal trend was observed for P max (Fig. 2C), which in this instance was associated with a significant interaction term (F 2,121 = 4.917, p = 0.009). There was a significantly steeper increase in P max with latitude under 10°C acclimation compared to 15°C and 20°C (Tukey Post-hoc comparison of regression slopes, p = 0.045 and 0.013, respectively) while there was no significant difference in slopes between 15 and 20°C acclimated snails (p = 0.933). Interestingly, P max and T peak were positively correlated across the entire dataset (T 126 = 3.555, p < 0.001, Pearson's correlation coefficient = 0.302), as were P max and CT max (T 126 = 4.765, p < 0.001, Pearson's correlation coefficient = 0.391), suggesting an association between  Fig. 1. Thermal performance curves. A) Graphical representation of the non-linear regression model used in our study, showing a TPC fitted to empirical heartrate data from an individual L. saxatilis (red circles). Curves are described using the functions t x b, with t describing the rising component of the TPC based on thermodynamic effects, and b representing the falling slope of the TPC due to protein instability. These functions incorporate four TPC parameters (Q 10 , C, T th , and d, see "Methods"), with T corresponding to ambient body temperature (°C). From these curves, five descriptors of thermal performance (P initial, P max , T breadth , T peak , and CT max ) are derived. See (Rezende and Bozinovic, 2019) for more information. (B\ \D) Averaged heart rates for each population during acute ramping assays. Heartrates were recorded at 1°C increments throughout ramping assays; error bars denote standard error. Curves were fitted to mean values using the "loess" function in the R package ggplot2. Note that the loess curves are presented as a visual representation of the raw data and are distinct from the non-linear regression approach shown in (A), which was used to fit TPCs for the main analysis. Visualisations of the non-linear model fits used in the main analysis are presented in supplementary fig. S2. maximum performance and upper limits of cardiac activity at the individual level.
Initial heart rates (P initial ) increased significantly with acclimation temperature (F 2, 122 = 5.352, p = 0.006) and across latitude (F 1, 122 = 7.463, p = 0.007), although there was no significant interaction (F 2, 122 = 2.451, p = 0.09; Fig. 2D). Post-hoc testing indicated that p Initial was highest in 20°Cacclimated individuals, followed by 15°Cacclimated individuals, and lowest under 10°Cacclimation (Tukey contrasts, all p < 0.05). A similar, albeit reversed pattern, was observed for the thermal sensitivity of heart rate, Q 10 , which decreased with increasing latitude (F 1, 122 = 5.540, p = 0.02) and increasing acclimation temperature (F 2, 122 = 4.011, p = 0.02; Fig. 2E). There was no significant interactive effect between latitude and acclimation temperature on Q 10 (F 2, 122 = 2.225, p = 0.112). Post-Hoc testing revealed significantly higher Q 10 values in 10°C acclimated snails compared to 15°C and 20°C (Tukey contrasts, both p < 0.001), but no significant difference between 15°C and 20°C acclimated snails (p = 0.999). In other words, while Q 10 was higher in warm-origin (low latitude) populations than in cold-origin populations, it was reduced under warm acclimation compared to cold acclimation, indicating contrasting effects of latitude vs. thermal acclimation on Q 10 . For T breadth , where latitude was treated as a categorical variable (see Methods), there was a significant interaction between latitude and acclimation temperature (F 6, 116 = 4.042, p = 0.01; Fig. 2F), although post-hoc testing revealed that this was chiefly driven by much higher T breadth in 15 and 20°Cacclimated Scottish snails, which differed significantly from all other populations at their respective acclimation temperatures (Tukey contrasts, all p < 0.05).

Principal component analysis
The first three components of our PCA encapsulated 87.8 % of total variation across performance metrics derived from the TPC fits (Q 10 , C, T th , logtransformed decay, T peak , CT max , T breadth and P max ). Plotting of PC1 against PC2 (39.7 % and 31.6 % of total variance respectively) reveals a high positive collinearity between the traits T peak , P max and CT max , as well as the descriptors T th and decay (Fig. 3), supporting the putative association between upper thermal limits and maximum performance suggested by our univariate analyses ( Fig. 2A-C). This cluster of variables was orthogonally separated from a second axis containing traits associated with thermal breadth, namely C, T breadth and the strongly negatively collinear Q 10 . PC3 (16.5 % of the remaining variation) was associated with variation in CT max and P max , which displayed a positive collinearity across our dataset (see Section 3.1). In general, TPCs for the different populations were largely displaced based on their thermal breadth, with the two higher latitude populations (Scotland and Iceland) displaying a shift towards higher thermal breadth and, in the case of the Icelandic population, also towards high T peak (grey arrows and accompanying annotations, Fig. 3). Acclimation to 10°C was generally associated with a shift towards lower thermal breadth, higher Q 10 , and reduced T peak , except in the Icelandic population, which instead displayed a pronounced shift towards high T peak relative to all other acclimation groups and populations.

Discussion
We assessed local adaptation and acclimation capacity of TPCs for cardiac performance in four latitudinally-separated intertidal populations of L. saxatilis from the Northeast Atlantic. Surprisingly, both the temperature at which cardiac performance was highest (T peak ) and the upper thermal limit (CT max ) increased with latitude, being higher in the northernmost population from Iceland and lowest in the southernmost population from Spain. Strikingly, populations from higher latitudes also exhibited elevated heartrates at T peak (P max ), as well as high basal heart rates (P initial ), when acclimated at low temperature (10°C). Because the high upper tolerance limits of northern populations run counter to the latitudinal gradient in temperature regimes experienced by the four populations, we posit that this pattern is unlikely to represent to thermal adaptation to cope with high temperatures. Instead, the close association between thermal limits and heart rate seen in our results suggests that elevated upper limits in high latitude populations are a consequence of maintaining a high metabolic machinery at low temperatures, a pattern associated with Metabolic Cold Adaptation (MCA). This is supported by the fact that warmacclimation (at 15 and especially 20°C) resulted in an overall reduction in maximum heart rates, especially in cold-adapted populations, potentially hinting at the existence of specific cold water enzyme pathways which are under selection at high latitudes. Our results suggest that metabolic performance at low temperatures and upper thermal limits are closely mechanistically linked, indicating that the interrelation between different aspects of the TPC may be a crucial factor in understanding variation in thermal stress responses across latitude.
A conceptual framework for interpreting our results is shown in Fig. 4. To compensate for lower thermodynamic reaction rates and reduced growing seasons in cold climates, cold-adapted populations are subject to thermal selection for elevated metabolic rates at low temperatures (Fig. 4A, greyed area) (Addo-Bediako et al., 2002;Leung et al., 2021;Villeneuve et al., 2021;White et al., 2012). This may, in turn drive other compensatory responses including higher growth and feeding rates (Conover et al., 2009). Due to kinetics, high metabolic rates at low temperatures translate into even higher metabolic rates under elevated temperatures (Gillooly et al., 2001). This elevated performance leads to an overall shift in the upper extremes of the TPC (correlated responses, Fig. 4A), as higher P max is positively associated with T peak and, ultimately, CT max (Angilletta, 2009;Angilletta et al., 2010;Frazier et al., 2006). Although high performance is costly to maintain for extended periods, it is likely that temperatures approaching T opt or CT max will rarely, if ever, be experienced in coldadapted populations, and may be further mitigated by the briefness of exposure as well as behavioural thermoregulation (Chapperon et al., 2017;Gunderson and Stillman, 2015;Miller and Denny, 2011). Thus, high performance in the upper regions of the TPC will not be actively selected against in these populations. By contrast, warm-adapted populations (Fig. 4B) must contend with higher average temperatures as well as with frequent and extended periods of heat stress during which metabolic costs rise dramatically. Thermal selection therefore acts to reduce the metabolic rate of warm-adapted organisms within their usual thermal range (Fig. 4B, greyed area). This has the benefit of reducing metabolic rates at high temperatures, including P max , thus making energy expenditure sustainable under prolonged heat exposure (Marshall and McQuaid, 2011), but conversely resulting in lower T peak and CT max due to thermodynamic constraints on the shape of TPCs.
Central to our conceptual model is the fact that local adaptation in the shape of TPCs presumably reflects a compromise between the competing effects of selection occurring at different temperatures, the relative importance of which will vary according to latitudinal differences in climatic regimes, as well as in the intensity of solar heating experienced during emersion ("non-climatic thermal stress"; Marshall et al., 2010). Although it does not factor in extreme microhabitat temperatures at low tide, our satellite temperature data (Fig. S1), nonetheless indicates substantial differences in the climatic regimes experienced by the four populations which may hint at selective pressures imposed upon the shape of TPCs at different latitudes. Notably, northern populations (especially the Icelandic population) frequently experience average temperatures of well below 10°C, meaning that responses in these populations, especially under cold acclimation, may be geared towards the maintenance of adequate metabolic rate at even lower temperatures. In such context, P initial and Q 10 may represent "adaptive" traits under strong selection, and indeed both traits displayed latitudinal shifts consistent with MCA in 10°C acclimated snails. Notably, Q 10 , was lower in cold-adapted populations, apparently contradicting the longstanding idea that populations from higher latitudes should display a greater thermal sensitivity of metabolism (Clarke, 1993;Hochachka and Somero, 2002). However, low Q 10 values may be advantageous in coldorigin populations based on first principles, as low Q 10 is associated with elevated metabolic rate at low temperature, as well as reducing the critical thermal minimum (Rezende and Bozinovic, 2019). By contrast, responses in the upper region of the TPC in northern populations (representing temperatures rarely experienced in nature) may have been unconstrained by thermal selection and therefore represented a "by-product" of selection at lower temperatures. For instance, if higher metabolic rate under MCA was associated with a switch to more flexible enzyme isoforms which remain stable across a wider range of temperatures, this may explain the higher P max , T opt and CT max , observed in the Icelandic population. Indeed, a previous study by Sokolova and Pörtner (2001) found evidence of higher respiratory enzyme activity and denaturation temperatures in a northern (White Sea) compared to a southern (North Sea) population of L. saxatilis.   Interestingly, under acclimation to 15 and 20°C, the Icelandic population used in our study displayed a shift in its response towards a pattern characteristic of the other populations studied, potentially representing a switch to enzyme isoforms also used by southern populations. This would indicate substantial metabolic plasticity in the Icelandic population, even when acclimated to temperatures outside the range they naturally experience on a regular basis. Latitudinal differences in thermal performance observed in this study could also be associated with morphological and ecological differences between the populations used, including genetic diversity resulting from their population history. Littorina saxatilis is known to exhibit high morphological diversity throughout its range (Reid, 1996), generally associated with microhabitat differences and including the formation of verticallysegregated ecotypes in some populations (Butlin et al., 2014). Although physiological differentiation has been associated with morphological variation across shore height (Dwane et al., 2021;Rolán-Alvarez et al., 1997;Sokolova et al., 2000), it is less clear how this may relate to differences across broader spatial scales such as latitude. In this study, all snails were sampled from upper shore populations belonging to the "moderate" ecotype sensu Reid (1996), characteristic of moderately wave-exposed shores. Despite this, morphological differences between the populations were present, notably the much larger body size in the Icelandic population. This may reflect differences in life history traits and energetics and potentially hints at other countergradient responses including faster growth rate in these populations, which may in turn be facilitated by the higher metabolic rate we observed in this population. Interestingly, Reid (1996) notes that large body size is characteristic of many high latitude populations of Littorina saxatilis (including northern Norway, the White Sea, Greenland and Canada), indicating that the response is not site-specific, although it is also suggested by Reid that this morphology may be due to a reduction in predation pressure due to the absence of Carcinus maenas. Unfortunately, little is known of growth trajectories and life-history in northern L. saxatilis populations, meaning it is unclear how higher metabolic rates may be linked to other aspects of the ecology of these populations, although this may provide a fruitful avenue for future research. Although we did not obtain genetic information from each of the four populations studied, phylogenetic studies (Doellman et al., 2011;Panova et al., 2011) suggest that the four populations also possess contrasting genetic histories which may have contributed to the physiological differences we identified. Notably, while the three northern populations of Littorina saxatilis used in the study migrated to their present-day locations since the end of the last glacial maximum (LGM), they did so from different ancestral populations, with those now inhabiting Iceland likely originating from glacial refugia in the Faroe Islands, while the English and Scottish populations likely originated in refugia located in the English channel or Irish sea. By contrast, the Spanish population was present on the NW Spanish coastline throughout the LGM (Panova et al., 2011). It is therefore possible that physiological differences we identified in this study may be the result of longstanding differences in the historical selective environments of the four populations throughout their evolutionary history, in addition to the effect of recent localised adaptation.
An additional factor potentially influencing differences in thermal selective regimes experienced by the populations used in our study is the degree of synchronicity between periods of emersion at low tide and exposure to midday summer temperatures (Finke et al., 2007;Helmuth et al., 2006;Lathlean et al., 2014;Sanford and Kelly, 2011). For example, Helmuth et al. (2006) found that Southern Californian sites may be protected from high temperatures due to low tide emersion generally occurring at night, while on a broader scale, certain locations (such as East Australia) exhibit particularly high occurrences of low tide aerial exposure during midday summers compared to comparable shorelines around the world (Finke et al., 2007). The extent of variation in synchronicity between tidal regimes and midday exposure is less well understood along the northeast Atlantic coastline, although in-situ temperature data has revealed that thermal regimes at certain European sites may differ from what would be anticipated based upon latitudinal location alone, perhaps also due to the influence of coastal upwellings (Seabra et al., 2011(Seabra et al., , 2015. A more extensive understanding of in-situ temperature regimes at northern sites, particularly in locations such as Iceland, would aid future work on thermal tolerance and metabolic cold adaptation in species with northerly distributions such as Littorina saxatilis. However, the clear latitudinal trend in average temperatures we observed from satellite data leads us to suspect that maximum temperatures experienced by the northern population used in this study would still be substantially lower than those experienced by the three southern populations. The close mechanistic association we observed between performance at low temperatures and upper thermal limits contrasts with patterns observed in tropical intertidal species, many of which display specializations to cope with extreme thermal stress at low tide (Marshall and McQuaid, 2011;Verberk et al., 2016;Wang et al., 2022). Under increasing temperatures, many such groups display a "thermally insensitive zone", a range of temperatures over which metabolic rate is held constant. This serves to Conceptual framework proposed to explain our results, demonstrating how correlated effects on the shape of TPCs may result from thermal selection on metabolism in A) cold-adapted and B) warm adapted populations. In both figures, the greyed area represents the range of temperatures over which thermal selection on metabolic rate is experienced by each population, while the black TPC represents a hypothetical ancestral population for the purposes of comparison. reduce energy exposure under prolonged thermal stress, resulting in a bimodal TPC characterized by an early peak associated with the maintenance of routine behaviour, and a much later peak presumably associated with the mobilization of the heat shock response prior to collapse (Verberk et al., 2016). In contrast to this study, bimodal peaks in such species appear to be mechanistically decoupled, suggesting that traits affecting performance at low temperatures vs. tolerance at high temperatures undergo independent thermal selection and are not interlinked (Marshall et al., 2010;Monaco et al., 2017). By contrast, the fact that cardiac activity in L. saxatilis increased steadily towards T peak in a manner consistent with our unimodal TPC model (Rezende and Bozinovic, 2019), may suggest that L. saxatilis lacks the specialist physiological adaptations present in tropical intertidal species to deal with extreme thermal stress, potentially because the range of temperatures over which metabolic depression occurs in tropical gastropodsfor example 33-45°C in Echinolittorina radiata (Wang et al., 2022)may lie outside the range of temperatures experienced by L. saxatilis for protracted periods. Behavioural thermal regulation, either through shell standing, aggregation, or reliance upon thermal refugia, is known to be prevalent in L. saxatilis (Chapperon et al., 2017;Dong et al., 2021Dong et al., , 2017 and the results of our study suggest that southern populations of L. saxatilis may rely upon behavioural mechanisms to a much greater extent than northern populations (Miller and Denny, 2011), potentially circumventing the need for physiological specialisms relied upon by tropical species. However, it is unclear to what extent reliance upon behavioural thermoregulation differs across the latitudinal range of L. saxatilis. Although the trait we refer to as T peak is often characterized as the thermal optimum (T opt ) in the TPC literature (Angilletta, 2009;Huey and Kingsolver, 1989), our findings highlight that the adaptive context implied by such terminology can be misleading when referring to TPCs for metabolism.
It is unclear to what extent our use of acute thermal ramping (6°C h −1 ) may have contributed to the close thermodynamic link we identified between performance at low temperatures and the upper regions of the TPC, as ramping rate can affect the limits of cardiac performance in intertidal species (Moyen et al., 2019). For instance, rapidly increasing temperatures may have led to a "snowball" effect on performance at higher temperatures, enhancing the correlated responses we saw in cold-adapted populations and potentially preventing metabolic depression from occurring in warm-adapted populations (Marshall and McQuaid, 2011;Verberk et al., 2016). It is possible that under slower rates of ramping, the Q 10 effects of MCA on upper regions of the TPC would have been reduced, leading to a diminishment of the countergradient latitudinal response seen in our results. However, it should be noted that 6°C h −1 represents an ecologically relevant ramping rate experienced during emersion for many temperate intertidal organisms (Bjelde and Todgham, 2013;Moyen et al., 2019), including L. saxatilis (C Dwane, pers. obs.), and indeed is slower than rates reported in some previous studies (Dowd et al., 2015;Logan et al., 2012;Moyen et al., 2019;Sokolova et al., 2000;Sokolova and Pörtner, 2003). In any case, given that the measurement of cardiac limits under acute ramping represent a widespread method in comparative physiological studies of intertidal organisms (Dong et al., 2021(Dong et al., , 2015Monaco et al., 2017), our study is significant in demonstrating that correlated responses in the shape of TPCs arising via thermal effects may play a role in influencing thermal sensitivity patterns across latitude, emphasising that it is important that future studies utilising approaches consider the potential impact of these effects. In addition, by highlighting the integrated nature of TPC traits, our study also provides a caution against the measurement of upper limits in isolation to infer local adaptation. We therefore advocate the use of approaches that enable variation across multiple TPC traits to be explored in tandem, both to better characterise patterns of intraspecific thermal adaptation across latitude and improve our ability to predict ecological responses to climate change at a global level.

Data availability
Data will be made available on request.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.