Supplementary Material

Heating tolerance of ectotherms is explained by temperature’s non-linear influence on biological rates

README
This document presents the full workflow of Heating tolerance of an ectotherm does not vary when temperature’s influence on biological rates is accounted for, including the full derivation of the equations, data processing and generating figures for the study. Figure and equation numbering here is independent of the main and supplemental text.

The analysis was conducted in R version 4.2.2 (2022-10-31 ucrt) (R Core Team 2022) running Tidyverse (Wickham 2019) & cowplot (Wilke 2020). This document was created using bookdown (Xie 2020) (Session info). Top right drop down menu can show code and download the Rmarkdown file.

---

Heating tolerance data

ARR <- readxl::read_excel("../data/Morley2018.xlsx",  skip = 8) %>% 
  # Add unique id
  rowid_to_column("row_ID") %>%  
  # remove blank rows
  filter(!is.na(Taxon)) %>%
  # rename columns
  rename("ta_high" = 11, #"High acclimation Temperature /°C" ,
         "ta_low" = 12, #"Low Acclimation Temperature/ °C" ,
         "ARR" = 17, #"CTMax ARR" ,
         "tc_high" = 16, #"CTMax/ °C", # CT max of the high acclimation temperature
         "acclim_window" = 13, #"Acclimation Window /°C",
         "heating_rate" = 22, #"°C change hr-1",
         "collection_latitude" = 19) %>%  #"Latitude of Collection/ °") %>%
  # fix typo in taxon
  mutate(Taxon = str_replace_all(Taxon, "amphibains", "amphibians")) %>%
  # Fix incorrect heating rate in DOI 10.1007/s10695-013-9793-7
  mutate(heating_rate = ifelse(row_ID  == 367, (0.3 * 60), heating_rate)) %>%
  # Calculate low temperature CT max
  mutate(tc_low = tc_high - (ARR * acclim_window)) 

We used the acclimation dataset available from the Polar Data Centre under an Open Government Licence v3.0 (Morley et al. 2018). The data was published in Morley et al. (2019). The dataset contains Acclimation Response Ratios (ARR), acclimation temperatures (\(T_a\)) and upper thermal tolerance limits (\(T_c\)) of ectotherms compiled from the literature between 1960 - 2015. Ectotherms were acclimated for between 1 and 300 days to a lower temperature (ta_low) or a higher temperature (ta_high) to give paired observations for each species. Animals were heated under a constant ramping temperature from \(T_a\) until the point of organismal collapse (\(T_c\)), called CT_max in Morley et al. (2019). Upper thermal tolerance limits were assayed for both acclimation temperatures in paired observations for a species (tc_low or tc_high). Within paired observations, each species was measured at the same constant heating rate, but heating rates varied among paired observations representing independent species or experimental assays. Data were only included if each species had \(T_c\) measured for two different starting \(T_a\), allowing us to compare heating tolerances both within and among species.

The dataset contained 369 observations of 316 species from 15 Classes representing marine, terrestrial and freshwater habitats (7 Phyla: Chordata, Arthropoda, Mollusca, Platyhelminthes, Echinodermata, Brachiopoda and Cnidaria).

The \(T_c\) of the lower acclimation temperature was not reported in Morley et al. (2019). We derived the \(T_c\) of the lower acclimation temperature from rearranging the formula for ARR. We first calculate the difference in \(T_c\) between the two acclimation temperatures from ARR multiplied by the difference in acclimation temperatures; called acclimation window (acclim_window) following Morley et al. (2019). We then calculate \(T_c\) of the lower acclimation temperature (tc_low) by subtracting the difference in \(T_c\) from the \(T_c\) of the higher acclimation temperature (tc_high).

\[ T_{c \space lower} = T_{c \space higher} - (ARR \times acclimation \space window) \]

---

Heating tolerance

We defined heating tolerance (\(H\)) as

\[\begin{equation} H = T_c - T_a \tag{1} \end{equation}\]

where \(T_a\) is acclimation temperature (°C, also starting temperature of the heating assay) and \(T_c\) is upper thermal tolerance limit (°C). We calculated heating tolerances (\(H\)) at the lower (ht_low) and higher acclimation (ht_high) temperature.

# Calculate heating tolerances, H (Equation 1)
ARR <- ARR %>% 
  mutate(ht_low = tc_low - ta_low,
         ht_high = tc_high - ta_high)

Figure 1: Heating tolerances at the lower and higher acclimation temperatures for paired observations on Log10 transformed axes.

---

Rescaling heating tolerance via heating duration

Physiological processes underlying heating responses scale non-linearly with temperature

# Define constants for rescaling
# Boltzmann's constant, k,  in units of eV K-1
k <- 8.61733 * 10 ^ -5
# Kelvin to Celsius conversion for normalisation T_K
Tk <- 273.15
# define \beta for a given activation energy, E = 0.6 (Equation 5)
beta_K <- 0.6 / (k * Tk)

If the endpoint of upper thermal tolerance limits is the product of biological processes and rates, then we may expect these rates to scale with temperature, independent of heating rate. Thus, acclimation temperature influences upper thermal tolerance limits via a non-linear thermal dependence. We use the Universal Temperature Dependence of biological rates (i.e., Boltzmann-Arrhenius equation, Equation (2)) to represent the non-linear effect of temperature on biological rates. We assume that the Boltzmann-Arrhenius equation adequately encompasses all physiological processes underlying heating tolerance and upper thermal tolerance limits, at least phenomenologically.

Mass-specific biological rates (\(r\)) increase with temperature following:

\[\begin{equation} r = R_0 e^{\frac{-E}{kT}} \tag{2} \end{equation}\]

Where \(R_0\) is a normalisation factor, \(E\) is activation energy, \(k\) is Boltzmann’s constant (\(8.617 {\times} 10^{-5}\) eV K^-1) and \(T\) is temperature in Kelvin. We use a value of 0.6 eV for \(E\) representing mean activation energy of biological rates across taxa (Dell, Pawar, and Savage 2011).

Although \(E\) varies greatly among taxa and biological rates of interest, varying the value of \(E\) does not meaningfully change our results (see Sensitivity of activation energy). \(R_0\) is not used in the calculation of rescaled heating duration and does not affect our results.

---

Time is implicit in heating tolerances through heating rates

Heating tolerances have a dimension of time (\(t\)) implicit in their definition. Under a constant heating rate, an organism must withstand heating over time between their acclimation temperature (\(T_a\)) and upper thermal tolerance limits (\(T_c\)).

In other words, heating tolerances also represent a duration of time (heating duration, \(\Delta t\)) between two time points (\(t_a\), the start of the heating assay, and \(t_c\), the end of the heating assay) that have a corresponding temperature (\(T_a\) and \(T_c\), respectively).

Therefore, as we can expect biological processes occurring during heating to scale with temperature we can also expect the durations of theses biological processes at one temperature, \(\Delta t(T)\), to scale with temperature. We can describe the relationship between heating durations starting from two temperatures (\(T_{a1}\) and \(T_{a2}\)) as:

\[\begin{equation} \frac{\Delta t(T_{a1})}{\Delta t(T_{a2})} \approx e^{\frac{E}{kT_{a1}}-\frac{E}{kT_{a2}}} \tag{3} \end{equation}\]

---

Rescaling the effect of acclimation temperature on heating durations

Following Equation (3), we can use a reference temperature to normalise heating durations for convenience and derive an expression that is independent of temperature but retains dimensions of time. We use 0°C as the reference temperature and define a normalisation constant (\(T_K\)) with a value of 273.15 (i.e. 0°C in Kelvin) for centring temperatures (\(T\)) in °C.

\[ \tau = \frac{T}{T_K} \]

Here, \(T\) are temperatures in degrees Celsius and \(\tau\) are normalised temperatures that are centred around the normalisation constant \(T_K\). \(\tau\) is unitless and can be converted back into degrees Celsius by \(T = \tau \times T_K\).

Substituting the normalisation constant \(T_K\) as \(T_{a1}\) and temperature \(T\) (in °C) as \(T_{a2}\) in Equation (3), we can derive an initial expression that rescales heating durations by accounting for acclimation temperature:

\[\Delta t_r = e^{\frac{E}{kT_K}}e^{-\frac{E}{kT}} \Delta t(T)\] Subscript \(r\) denotes variables rescaled by the non-linear temperature dependence of physiological rates (Table 1). \(\Delta t_r\) are rescaled heating durations that correspond to non-rescaled heating duration (\(\Delta t\)).

Rescaled heating duration can be written as an integral between the start (\(t_a\)) and end times (\(t_c\)) of the heating assay, assuming a constant temperature ramp:

\[\begin{equation} \Delta t_r = e^{\frac{E}{kT_K}} \int_{t_a}^{t_c} e^{-\frac{E}{kT(t)}} dt \tag{4} \end{equation}\]

We can simplify Equation (4) by using the normalised temperature scale (\(\tau\)) and removing the dimension of temperature.

On a normalised temperature scale (\(\tau\)), \(\frac{E}{kT}\) in Equation (4) becomes:

\[\begin{align} \frac{E}{kT} &= {\frac{E}{kT_K}} {\frac{1}{1 + \frac{T}{T_K}}} \\ & = {\frac{E}{kT_K}} {\frac{1}{1 + \tau}} \\ & \approx \beta (1-\tau) \end{align}\]

where \[\begin{equation} \beta = \frac{E}{k T_K} \tag{5} \end{equation}\]

\(\beta\) is a constant with an approximate value of 25.5 using \(E =\) 0.6 (beta_K).

Thus Equation (4) simplifies to

\[\begin{equation} \Delta t_r = \int_{t_a}^{t_c} e^{\beta \tau (t)} dt \tag{6} \end{equation}\]

describing rescaled heating duration as a function of time. Rescaled heating duration (\(\Delta t_r\)) can be considered as the effective elapsed duration of the heating assay once the non-linear effect of acclimation temperature has been accounted for.

Equation (6) can equivalently be expressed as an integration over normalised temperature scales (\(\tau\)).

\[\begin{equation} \Delta t_r = \frac{1}{\lambda} \int_{\tau_{a}}^{\tau_{c}} e^{\beta \tau} d \tau \tag{7} \end{equation}\]

Here, \(\tau_a\) is acclimation temperature and \(\tau_c\) is upper thermal tolerance limit on normalised temperature scales and correspond with \(T_a\) and \(T_c\) on degree Celsius scales, respectively. \(\lambda\) is heating rate (°C) on a normalised temperature scale, where \(\lambda = Heating \space rate / T_K\) with units of time^-1.

---

Calculating rescaled heating duration from heating tolerance

Integration of Equation (7) gives an expression for rescaled heating duration \(\Delta t_r\) that describes the effect of the non-linear change in physiological rates as temperatures increased during the heating assay, on normalised temperature scales (Equation (8)).

\[\begin{equation} \Delta t_r = g_\beta (\Delta \tau) \frac{e^{\beta \tau_a}}{\lambda} \tag{8} \end{equation}\]

\(\Delta \tau\) is heating tolerance on the normalised temperature scale, where \(\Delta \tau = \tau_c - \tau_a\). \(\Delta \tau\) is not corrected for the non-linear temperature dependence of physiological rates. The function \(g_\beta\) accounts for the change in rescaled temperature in the heating assay via integration:

\[\begin{align} g_\beta (x) &= \int_{\tau_{a}}^{\tau_{c}} e^{\beta x} dx \\ &= \frac{e^{\beta x} - 1}{\beta} \tag{9} \end{align}\]

Thus, Equation (8) is used to calculate rescaled heating duration \(\Delta t_r\) from a known heating tolerance (as \(H\) or \(\Delta \tau\)). Normalised heating tolerances \(\Delta \tau\) can readily be converted into degrees Celsius to calculate \(H\) or vice versa. Rescaled heating durations \(\Delta t_r\) measured at two acclimation temperatures are comparable within a species.

---

Calculating heating tolerance from rescaled heating duration

Conversely, normalised heating tolerance \(\Delta \tau\) can be calculated from rescaled heating duration \(\Delta t_r\) using the inverse of Equation (9).

\[\begin{equation} \Delta \tau = f_\beta (e^{-\beta \tau_{a}} \lambda \Delta t_r) \tag{10} \end{equation}\]

Where \(f_\beta\) is the inverse of \(g_\beta\) (Equation (9)) that also accounts for the change in rescaled temperature in the heating assay.

\[\begin{align} f_\beta(x) &= \frac{\log (1 + \beta x)}{\beta} = : g^{-1}_\beta(x) \tag{11} \end{align}\]

Equations (9) and (11) approach the identity for \(\beta x \ll 1\).

Thus, Equation (10) can be used to calculate normalised heating tolerances \(\Delta \tau\) from a known rescaled heating duration \(\Delta t_r\), and thus can be used to calculate heating tolerance in degrees Celsius, \(H\).

Equations (8) and (10) involves two steps: first the heating rate is rescaled by the term \(e^{\beta \tau_{a}}\) that depends on the starting temperature of the assay (to correct for acclimation temperature), second Equations (9) and (11) account for the change in temperature as the heating assay progresses via integration of the Universal Temperature Dependence.

Rescaling heating durations and heating tolerance is dependent on three assumptions: 1) Heating tolerance is arbitrarily defined with respect to 0°C, 2) heating rate is constant, and 3) the Arrhenius relationship is not affected by the physiological state of the organism or the experimental conditions (i.e., \(E\) does not change).

# function for delta t_r that uses g_beta_fun (Equation 8)
delta_tr_fun <- function(tau_a, tau_c, lambda){
  exp(beta_K * tau_a) / lambda * g_beta_fun(tau_c - tau_a)
}

# function g of x and the constant beta (Equation 9)
g_beta_fun <- function(x){(exp(beta_K * x) - 1) / beta_K}

# function f of x and the constant beta (Equation 11)
f_beta_fun <- function(x){log(1 + beta_K * x) / beta_K}

We calculated heating tolerances at the two acclimation temperatures on the normalised temperature scale. We then apply Equation (8) to calculate rescaled heating duration \(\Delta t_r\).

# normalise Ta, Tc and heating rate by Tk, calculate normalised heating tolerance
ARR <- ARR %>% 
  mutate(tau_a_low = ta_low / Tk,
         tau_a_high = ta_high / Tk,
         tau_c_low = tc_low / Tk,
         tau_c_high = tc_high / Tk,
         lambda = heating_rate / Tk,
         delta_tau_low = tau_c_low - tau_a_low,
         delta_tau_high = tau_c_high - tau_a_high)

# Calculate rescaled heating duration for lower acclimation temperatures (Equation 8)
ARR <- ARR %>% 
  mutate(delta_tr_low = pmap_dbl(list(tau_a = tau_a_low,
                                  tau_c = tau_c_low,
                                  lambda = lambda),
                             delta_tr_fun))

# Calculate rescaled heating duration for higher acclimation temperatures (Equation 8)
ARR <- ARR %>% 
  mutate(delta_tr_high = pmap_dbl(list(tau_a = tau_a_high,
                                   tau_c = tau_c_high,
                                   lambda = lambda),
                              delta_tr_fun))

Rescaled heating durations \(\Delta t_r\) are influenced by heating rate with larger rescaled heating durations associated with slower heating rates (Figure 2). This dependency with heating rate is explored below.

Figure 2: Rescaled heating duration at the lower and higher acclimation temperature on Log10 axes have a strong dependency with heating rate (colours).

---

Dimensionless rescaled heating duration

A simple way of visualising rescaled heating duration without its dependency on heating rate is to multiply rescaled heating duration by heating rate to give a dimensionless quantity (\(H_r\)). This dimensionless quantity demonstrates the strong relationship between rescaled heating durations at either acclimation temperature, regardless of heating rate (Figure 3). However, this dimensionless analysis has limited utility in understanding temperature-rate-time dynamics compared with \(\Delta t_r\). \(\Delta t_r\) includes the integration of the non-linear temperature dependency of biological rates and is our main aim of our analysis.

Figure 3: Dimensionless rescaled heating durations at the lower and higher acclimation temperatures on Log10 axes.

---

Sensitivity of activation energy

Intraspecific activation energy of various thermally-dependent traits ranges between 0.2 and 1.2 with a median of 0.55 eV (Dell, Pawar, and Savage 2011). Changing activation energy and thus \(\beta\) does not meaningfully change the reported patterns (Figure 4).

delta_tr_gen is a generalised function to rescale heating tolerance with \(E\) as a variable. \(\beta\) values are calculated for each value of \(E\).

# Generalised function
delta_tr_gen <- function(x, ev, k = 8.61733*10^-5, Tk = 273.15) { # ev is activation energy
  
  # calculate beta K constant
  beta_K_sensi <- ev / (k * Tk)
  
  # function for delta tr that uses g_beta_fun (Equation 8)
  delta_tr_gen <- function(tau_a, tau_c, lambda){
    exp(beta_K_sensi * tau_a) / lambda * g_beta_fun_gen(tau_c - tau_a)
  }

  # function f of x and the constant beta
  f_beta_fun_gen <- function(x){log(1 + beta_K_sensi * x) / beta_K_sensi}
  
  # function g of x and the constant beta
  g_beta_fun_gen <- function(x){(exp(beta_K_sensi * x) - 1) / beta_K_sensi}
  
  # normalise Ta, Tc and heating rate by Tk
  x <- x %>% # Calculate heating rate as lambda
    mutate(tau_a_low = ta_low / Tk,
           tau_a_high = ta_high / Tk,
           tau_c_low = tc_low / Tk,
           tau_c_high = tc_high / Tk,
           lambda = heating_rate / Tk,
           delta_tau_low = tau_c_low - tau_a_low, # calculate delta taus
           delta_tau_high = tau_c_high - tau_a_high)

  # map over low temperature experiments
  x <- x %>% 
    mutate(delta_tr_low = pmap_dbl(list(tau_a = tau_a_low,
                                    tau_c = tau_c_low,
                                    lambda = lambda),
                               delta_tr_gen))
  
  # do the same with the high temp acclimated data
  x <- x %>% 
    mutate(delta_tr_high = pmap_dbl(list(tau_a = tau_a_high,
                                     tau_c = tau_c_high,
                                     lambda = lambda),
                                delta_tr_gen))
  
  # Add predictions of high from low and conversely low from high
  x <- x %>% 
    mutate(activation = ev,
           delta_tau_high_hat = f_beta_fun_gen(exp(beta_K_sensi * (tau_a_low - tau_a_high)) * g_beta_fun_gen(delta_tau_low)), # predict high from low
           delta_tau_low_hat = f_beta_fun_gen(exp(beta_K_sensi * (tau_a_high - tau_a_low)) * g_beta_fun_gen(delta_tau_high))) # predict low from high

  return(x) # as list
}

The dataset is replicated for a vector of activation energies, then the delta_tr_gen function is applied to each element of the list using mapply.

# define activation energies
activation <- seq(0.2, 1.2, 0.2)
# repeat ARR dataset over length of activation energies
ARR_list <- replicate(length(activation), ARR, simplify = FALSE)
# integrate under curve, merge to single data frame
ARR_list <- do.call(rbind, mapply(ARR_list, FUN = delta_tr_gen, ev = activation, SIMPLIFY = FALSE))

# plot
sensi_plot <- ggplot(data = ARR_list, 
                 mapping = aes(x = delta_tr_low, 
                               y = delta_tr_high)) + 
  geom_point(col = "grey") + 
  geom_smooth(method = "lm", se = FALSE, col = "black") + 
  geom_abline(slope = 1, intercept = 0, lty = 2, alpha = 0.6) +
  scale_x_continuous(expand = c(0.01,0), trans = "log10", labels = function(x) format(x, scientific = FALSE)) +
  scale_y_continuous(expand = c(0.01,0), trans = "log10", labels = function(x) format(x, scientific = FALSE)) +
  # scale_colour_viridis_c(name = "Activation energy") +
  facet_wrap(.~activation, ncol=3, scales = "free") +
  labs(x = expression("Lower temperature rescaled heating duration"), y = expression("Higher temperature rescaled heating duration")) +
  theme_classic(base_size = 12)

sensi_plot

Figure 4: Variation in activation energy (eV). Observations (points) with a fitted linear model (solid line) fall along the identity line (dashed line).

---

Predicting heating tolerance for a new acclimation temperature

We can use rescaled heating durations to predict heating tolerances measured at a different acclimation temperature (\(a2\)) from a known acclimation temperature (\(a1\)) and heating tolerance.

# Add predictions of high from low and conversely low from high
ARR <- ARR %>% 
  mutate(delta_tau_high_hat = f_beta_fun(exp(beta_K * (tau_a_low - tau_a_high)) * g_beta_fun(delta_tau_low)), # predict high from low (Equation 12)
         delta_tau_low_hat = f_beta_fun(exp(beta_K * (tau_a_high - tau_a_low)) * g_beta_fun(delta_tau_high))) # predict low from high (Equation 13)

Assuming heating tolerances are the same

The simplest model to predict a new heating tolerance from a new acclimation temperature is to assume that heating tolerances at the higher and lower acclimation temperature experiments are the same. For example, we can predict the heating tolerance at the higher acclimation temperature \(\Delta \tau (\tau_{a2})\), from the lower acclimation temperature, \(\Delta \tau (\tau_{a1})\).

\[\begin{equation} \Delta \tau (\tau_{a2}) \approx \Delta \tau (\tau_{a1}) \tag{12} \end{equation}\]

pred1 <- ggplot(data = ARR,
       mapping = aes(x = delta_tau_low * Tk,
                     y = delta_tau_high * Tk)) + 
  geom_point() + 
  geom_abline(intercept = 0, slope = 1, lty = 2, alpha = 0.6) + 
  geom_smooth(method = "lm", se = FALSE, col = "darkgrey") + 
  labs(x = expression("Predicted heating tolerance (\u00B0C)"),
       y = expression("Observed heating tolerance (\u00B0C)")) + 
  scale_x_continuous(expand = c(0,0), limits = c(0, 40)) +
  scale_y_continuous(expand = c(0,0), limits = c(0, 40))

pred1

Figure 5: Predicting higher acclimation temperature heating tolerance from lower acclimation temperature heating tolerance assuming heating tolerance is the same.

Assuming heating tolerance is the same at both acclimation temperatures overestimates the heating tolerance at the higher temperature (Figure 5).

Correct for acclimation temperature difference

Next, we can correct for the difference in acclimation temperatures by rescaling by the acclimation temperature via \(e^{\beta (\tau_{a1} - \tau_{a2})}\):

\[\begin{equation} \Delta \tau (\tau_{a2}) \approx \Delta \tau (\tau_{a1}) e^{\beta (\tau_{a1} - \tau_{a2})} \tag{13} \end{equation}\]

pred2 <- pred1 + aes(x = (delta_tau_low * exp(beta_K * (tau_a_low - tau_a_high))) * Tk) +
  labs(x = expression("Predicted heating tolerance (\u00B0C)")) 
pred2

Figure 6: Predicting higher acclimation temperature heating tolerance correcting for acclimation temperature

Using Equation (13) underestimates the heating tolerance at the higher acclimation temperature (Figure 6).

Account for difference in acclimation temperature and change in temperature

Finally, we can derive a more precise prediction that fully accounts for the change to heating tolerance that occurs during the experiment as the temperature is increased using Equation (14).

\[\begin{equation} \Delta \tau (\tau_{a2}) \approx f_\beta \big(e^{\beta (\tau_{a1} - \tau_{a2})} g_\beta(\Delta \tau(\tau_{a1})) \big) \tag{14} \end{equation}\]

pred3 <- pred1 + aes(x = delta_tau_high_hat * Tk) + 
  labs(x = expression("Predicted heating tolerance (\u00B0C)"))
pred3

Figure 7: Correction for acclimation temperature and change in temperature over time to predict heating tolerance at the higher acclimation temperature from the lower acclimation temperature.

The \(R^2\) is 93.45%. And as we can predict heating tolerance at the higher acclimation temperature from the lower acclimation temperature in Figure 7, so we can also predict heating tolerance at the lower acclimation temperature from heating tolerance at the higher temperature using Equation (15) (Figure 8).

\[\begin{equation} \Delta \tau (\tau_{a1}) \approx f_\beta \big(e^{\beta (\tau_{a2} - \tau_{a1})} g_\beta(\Delta \tau(\tau_{a2})) \big) \tag{15} \end{equation}\]

pred4 <- pred1 + aes(x = delta_tau_low_hat * Tk,
                     y = delta_tau_low * Tk) + 
  # ggtitle("Predict low from high") +
  labs(x = expression("Predicted heating tolerance (\u00B0C)"), 
       y = expression("Observed heating tolerance (\u00B0C)")) 

pred4

Figure 8: Predicting lower temperature heating tolerance from higher temperature heating tolerance

The \(R^2\) is 91.24%. Equations (14) and (15) converts normalised heating tolerance at one acclimation temperature into rescaled heating tolerance (via \(g_\beta\)), corrects for the difference in acclimation temperature along the non-linear Universal Temperature Dependence (via \(e^{\beta\tau}\)), and converts rescaled and biological rate-corrected heating tolerance into normalised heating tolerance at the other acclimation temperature (via \(f_\beta\)).

---

The effect of heating rate on rescaled heating duration

varying_rates <- read.csv("Morley_et_al_2016.csv") %>% 
  rename(heating_rate = 8,
         ta = 9,
         tc = 10) %>% 
   mutate(Genus = ifelse(Species == "borchgrevinki", "Pagothenia", genus)) %>% # Synonym Trematomus Pagothenia Pagothenia in ARR dataset
  mutate(Name = paste(Genus, Species), 
         experiment = "variable",
         lambda = heating_rate / Tk) %>%
  mutate(Phyla = case_when(phyla == "mollusc" ~ "Mollusca",
                           phyla == "Fish" ~ "Chordata",
                           phyla == "Ascidians" ~ "Chordata",
                           phyla == "Crustacean" ~ "Arthropoda",
                           TRUE ~ phyla)) %>%
  group_by(Name) %>%
  distinct() %>% # remove duplicated rows
  filter(n()>1) %>% # remove species with 1 observation
  ungroup()

To examine the relationship between heating rate and rescaled heating duration further, we used a second dataset containing \(T_c\) measured from the same \(T_a\) under different heating rates (\(\lambda\)) (Morley et al. 2016). Species with only one observation are removed (n = 3). The dataset has 107 observations from 7 phyla. Trematomus borchgrevinki is renamed to Pagothenia borchgrevinki because Pagothenia is the accepted genus and is also used in the ARR dataset. We combined this dataset with data from the constant heating rate experiments for species that also had variable heating rate data with them (129 observations from 37 species).

The original ARR dataset was converted to a long format to join with the varying heating rate dataset (species_data). The datasets are identified by experiment (variable or constant) and occurrence (identifying species that occur in both datasets).

species_data <- ARR %>%
  # Transform the ARR dataset to long format
  pivot_longer(cols = starts_with(c("ta_", "tc_")),
               names_to = c(".value", "assay"),
               names_sep = "_") %>%
  mutate(experiment = "constant") %>% 
  # Join datasets and get occurences of species in each dataset
  full_join(.,
            varying_rates, 
            by = c("Genus", "Species", "Name", "experiment", "ta", "tc", "lambda"),
            suffix = c("_constant", "_variable")) %>% 
  # Normalise temperature
  mutate(tau_a = ta / Tk, 
         tau_c = tc / Tk,
         delta_tau = tau_c - tau_a) %>% 
  group_by(Name) %>% 
  mutate(occurrence = length(unique(experiment))) %>%
  ungroup() %>% 
  # Keep species with more than 1 heating rate from either dataset
  filter(experiment == "variable" | occurrence == 2) %>% 
  # Calculate rescaled heating duration (Equation 8)
  mutate(delta_tr = pmap_dbl(list(tau_a = tau_a,
                                  tau_c = tau_c,
                                  lambda = lambda),
                             delta_tr_fun))

rate_ratios <- species_data %>% 
  # Group by name
  group_by(Name) %>% 
  # For each species create all combinations of heating rate with covariates
  group_modify(~expand_grid(.x, 
                            select(.x, lambda, delta_tr, delta_tau, tau_a),
  # Add numbered suffix to all columns; unique for duplications
                            .name_repair = function (x) ave(x, x, FUN = function(i) str_c(i, seq_along(i))))) %>%
  # Calculate ratio of rescaled duration and heating rate (Equation 19)
  mutate(lambda_ratio = lambda2/lambda1,
         t_ratio = delta_tr2/delta_tr1) %>%
  # Calculate scaling coefficient gamma for each species (Equation 19)
  group_modify(function(x, y){
    x$gg <- abs(as.numeric(coef(lm(log(t_ratio) ~ log(lambda_ratio), x))[2]))
    return(x)}) %>% 
  distinct() %>% # Remove duplications created by expand grid
  ungroup() %>% 
  # Calculate predicted normalised heating tolerance at the second acclimation temperature, corrected for heating rate (Equation 20)
  mutate(delta_tau2_hat = f_beta_fun(exp(beta_K * (tau_a1-tau_a2)) * (lambda_ratio ^ (1 - gg)) * g_beta_fun(delta_tau1)),
         delta_tau2_hat_inter = f_beta_fun(exp(beta_K * (tau_a1-tau_a2)) * (lambda_ratio ^ (1 - 0.7656548)) * g_beta_fun(delta_tau1)))

# interspecific gamma
gamma <- abs(coef(glm(log(t_ratio) ~ log(lambda_ratio), data = rate_ratios))[2])

---

Testing a null effect of heating rates on rescaled heating duration

We first considered a null hypothesis that rescaled heating duration does not depend on heating rate. Equation (10) suggests that for any heating rate:

\[\begin{equation} \Delta \tau (\lambda)= f_\beta (e^{-\beta \tau_{a}} \lambda \Delta t_r(\lambda)) \tag{16} \end{equation}\]

If the null hypothesis is supported, then \(\Delta t_r (\lambda) = \Delta t_r (0)\) and thus Equation (16) becomes

\[\begin{equation} \Delta \tau (\lambda) = f_\beta (e^{-\beta \tau_{a}} \lambda \Delta t_r (0)) \tag{17} \end{equation}\]

which is a logarithmic, temperature-dependent relationship between observed heating tolerance and heating rate. Equation (17) can be tested by conducting ramping assays starting at the same acclimation temperature but using different heating rates. From these data we can infer \(\Delta t_r (0)\) from the experiment with the slowest heating rate \(\lambda_{min}\). We should then have that

\[\begin{equation} \Delta \tau (\lambda) = f_\beta (\hat{\lambda}), \ \text{where} \ \hat{\lambda} = \frac{\lambda}{\lambda_{min}} g_\beta (\Delta \tau (\lambda_{min})) \tag{18} \end{equation}\]

In a plot of normalised heating tolerance (\(\Delta \tau\)) against the transformed variable \(f_\beta (\hat{\lambda})\), the null hypothesis is supported if the points fall along the identity (1:1) line. If points systematically fall below the 1:1 line, then the null hypothesis is not supported.

g_departure <- species_data %>%
  filter(experiment == "variable") %>% 
  group_by(Name) %>% 
  mutate(lambda_min = min(lambda),
         delta_tau_min = delta_tau[which.min(lambda)]) %>% 
  ungroup() %>% 
  ggplot(mapping = aes(x = f_beta_fun(lambda / lambda_min * g_beta_fun(delta_tau_min)) * Tk,
                       y = delta_tau * Tk,
                       col = Name)) + 
  geom_point() + 
  geom_line(alpha = 0.5) +
  geom_abline(slope = 1, intercept = 0, lty = 2, alpha = 0.6) +
  labs(x = expression("Predicted heating tolerance (\u00B0C)"), 
       y = expression("Observed heating tolerance (\u00B0C)")) +
  scale_colour_viridis_d(option = "C") +
  scale_x_continuous(expand = c(0,0), limits = c(0, 100)) +
  scale_y_continuous(expand = c(0,0), limits = c(0, 25)) +
  theme(legend.position="none", plot.margin = margin(5.5, 7, 5.5, 5.5, unit = "pt"))

g_departure

Figure 9: Predicted heating tolerance versus observed heating tolerance for species with varying heating rates (107 observations). Colours are species (n = 37).

We used the slowest heating rate of each species to test whether heating rate had no effect on rescaled heating duration and thus predicted heating tolerance. The null hypothesis was not supported because points systematically fall below the 1:1 line (Figure 9).

To visualize the effect of heating rate on rescaled heating tolerance, we plot the ratio of rescaled heating duration for a paired observation within a species with the corresponding ratio for paired heating rate. The relationship between the ratio of rescaled heating duration and the ratio of heating rate is described by an allometric scaling equation in the form \(\Delta t_r (\lambda) \sim c \lambda^{- \gamma}\) where \(c\) depends only on the experimental protocol and the species-specific scaling exponent, \(\gamma\). Points should fall along the identity line if heating rate does not affect rescaled heating tolerance (i.e., a scaling exponent of -1).

Change in rescaled heating duration, and thus heating tolerance, scales with change in heating rate according to an allometric scaling equation. There is some species-specific variation in the slope of the regression \(\gamma\) (colours). Overall, there is an interspecific scaling exponent of approximately \(\gamma =\) 0.766 calculated from a linear regression of ratios of heating rate and ratios of rescaled heating duration on Log₁₀ scales (Figure 10).

g_scaling <- ggplot(data = rate_ratios,
                    mapping = aes(x = lambda_ratio, 
                                  y = t_ratio, 
                                  group = Name,
                                  col = Name)) + 
  geom_line(stat ="smooth", method = "lm", alpha = 0.5) +
  geom_point() + 
  geom_abline(slope = -1, intercept = 0, lty = 2, alpha = 0.6) +
  geom_smooth(aes(group = 1), method = "lm", se = FALSE, col = 1) +
 scale_x_log10(
   breaks = scales::trans_breaks("log10", function(x) 10^x),
   labels = scales::trans_format("log10", scales::math_format(10^.x))
 ) +
 scale_y_log10(
   breaks = scales::trans_breaks("log10", function(x) 10^x),
   labels = scales::trans_format("log10", scales::math_format(10^.x))
 ) +
  ggtitle(paste0("slope = ", 
                 round(gamma, 3))) +
  labs(x = "Heating rate ratio", y = "Rescaled heating duration ratio") +
  scale_colour_viridis_d(option = "C") +
  theme(legend.position="none")

g_scaling

Figure 10: Interspecific scaling of rescaled heating duration with heating rate. Variation within each species is shown by the fitted lines between points (129 observations). The dashed identity line is the expected relationship following the null hypothesis. Axes are Log10 transformed. Colours are species (n = 37).

---

Quantifying the effect of heating rates on rescaled heating tolerance

Since higher heating rates reduces rescaled heating tolerance following a phenomenological power law, then the relationship between two experiments done at two heating rates, \(\lambda_1\) and \(\lambda_2\), regardless of the starting temperature, can be described by:

\[\begin{equation} \frac{\lambda_2 \Delta t_r (\lambda_2)}{\lambda_1 \Delta t_r (\lambda_1)} \approx \bigg(\frac{\lambda_2}{\lambda_1}\bigg)^{1-\gamma} \tag{19} \end{equation}\]

and therefore that from one heating tolerance \(\Delta \tau_1\) we should be able to predict any other \(\Delta \tau_2\). Using Equation (19) we can deduce \(\Delta t_r (\lambda_2)\) and from Equations (8) and (10) leads to:

\[\begin{equation} \Delta\tau_2 = f_{\beta}\Bigg(e^{\beta(\tau_{a1}-\tau_{a2})}\bigg(\frac{\lambda_2}{\lambda_1}\bigg)^{1-\gamma}g_{\beta}\big(\Delta\tau_1\big)\Bigg) \tag{20} \end{equation}\]

That is, we can predict one heating tolerance \(\Delta\tau_2\) from another \(\Delta\tau_1\) using Equation (20) which incorporates the correction for the effect of heating rate on rescaled heating duration using a species specific scaling exponent (\(\gamma\)). We use species-specific scaling exponents in subsequent calculations for precision. Doing so, we see strong concordance of the observed heating tolerance (delta_tau2) and the estimated value using Equation (20) (delta_tau2_hat).

g_predict_tau <- rate_ratios %>% 
  filter(t_ratio != 1 & lambda_ratio != 1) %>% # remove trivial pairs where it is the lowest heating rate
  ggplot(mapping = aes(x = delta_tau2_hat * Tk, 
                       y = delta_tau2 * Tk,
                       col = Name)) + 
  geom_abline(slope = 1, intercept = 0, lty = 2, alpha = 0.6) +
  geom_point() + 
  geom_smooth(method = "lm", se = FALSE, col = 1) + 
  labs(x = expression("Predicted heating tolerance (\u00B0C)"), 
       y = expression("Observed heating tolerance (\u00B0C)")) +
  scale_colour_viridis_d(option = "C") +
  scale_x_continuous(expand = c(0,0), limits = c(0, 30)) +
  scale_y_continuous(expand = c(0,0), limits = c(0, 30)) +
  theme(legend.position="none")

g_predict_tau

Figure 11: Predicted heating tolerance corrected using species-specific scaling exponents for the allometric relationship between rescaled heating duration and heating rate. Colours are species.

We can plot the same as Figure 11 but using the interspecific value of 0.766 for \(\gamma\) (Figure 12). Doing so increases our prediction error for a species.

Figure 12: Predicted heating tolerances calculated using a value of 0.766 for gamma, representing an interspecific parameter. Colours are species.

---

Generalised graphs of rescaled heating tolerances

We can visualise the general relationship described by Equation (20), using median values of heating rate, \(T_c\) and \(T_a\) across all species and \(\gamma =\) 0.766.

sp_heat_rate = median(ARR$heating_rate)
sp_ta = median(ARR$ta_low)
sp_tc = median(ARR$tc_low)

heating_rate <- round(sort(c(10^(seq(log10(0.004), log10(60), length.out = 14)), 0.015, 1, 0.042, 18)), 3) # on log10 scale for equidistant points
heating_rate <- heating_rate[!heating_rate %in% c(0.018, 0.037, 13.667)] # neaten lines
ta <- seq(5, 40, 1) # range of acclimation temperatures
pal_cols <- viridisLite::viridis(n = length(heating_rate), end = 1) # colour palette for points

gen_ecto <- expand_grid(heating_rate, ta) %>% 
  mutate(sp_lambda = sp_heat_rate / Tk,
         sp_tau_a = sp_ta / Tk,
         sp_tau_c = sp_tc / Tk,
         sp_delta_tau = sp_tau_c - sp_tau_a,
         tau_a = ta / Tk, 
         delta_tau_a =  sp_tau_a - tau_a,
         lambda = heating_rate / Tk,
         lambda_ratio = lambda/sp_lambda)  %>% 
  mutate(delta_tau_hat =  f_beta_fun(exp(beta_K * delta_tau_a) * (lambda_ratio ^ (1 - gamma)) * g_beta_fun(sp_delta_tau))) %>% 
  ggplot() + 
  geom_line(aes(ta, delta_tau_hat * Tk, group = heating_rate, col = heating_rate)) +
  labs(x = expression("Acclimation temperature (\u00B0C)"), 
       y = expression("Expected heating tolerance (\u00B0C)")) +
  scale_colour_viridis_c(name = expression("Heating rate (\u00B0C h" ^-1 * ")"), trans = "log10", labels = function(x) signif(x)) 

gen_ecto

Figure 13: Generalised graph of predicted heating tolerance by acclimation temperature and heating rate (colours).

Figure 13 can also be visualised as predicted upper thermal tolerance limit \(T_c\) for any acclimation temperature \(T_a\).

expand_grid(heating_rate, ta) %>% 
  mutate(sp_lambda = sp_heat_rate / Tk,
         sp_tau_a = sp_ta / Tk,
         sp_tau_c = sp_tc / Tk,
         sp_delta_tau = sp_tau_c - sp_tau_a,
         tau_a = ta / Tk, 
         delta_tau_a =  sp_tau_a - tau_a,
         lambda = heating_rate / Tk,
         lambda_ratio = lambda/sp_lambda)  %>% 
  mutate(delta_tau_hat =  f_beta_fun(exp(beta_K * delta_tau_a) * (lambda_ratio ^ (1 - gamma)) * g_beta_fun(sp_delta_tau))) %>% 
  ggplot() + 
  geom_line(aes(ta, delta_tau_hat * Tk + ta, group = heating_rate, col = heating_rate)) +
  labs(x = expression("Acclimation temperature (\u00B0C)"), 
       y = expression("Upper thermal tolerance limit (\u00B0C)")) +
  scale_colour_viridis_c(name = expression("Heating rate (\u00B0C h" ^-1 * ")"), trans = "log10", labels = function(x) signif(x)) 

Figure 14: Generalised graph of predicted upper thermal tolerance limit by acclimation temperature and heating rate (colours).

---

Notation table

Table 1: Main mathematical notation and description used.

Symbol	Description
\(H\)	Heating tolerance; range of temperatures between \(T_a\) and \(T_c\)
\(\Delta t\)	Elapsed duration of a heating assay between \(t_a\) and \(t_c\)
\(_r\)	Variables rescaled by the Universal Temperature Dependence
\(_a\)	Variables referring to acclimation or starting temperature of a heating assay
\(_c\)	Variables referring to the upper thermal tolerance limit of a heating assay
\(\tau\)	Normalised temperature scale centred around \(T_K\) (normalisation constant 273.15K)
\(\lambda\)	Normalised heating rate centred around \(T_K\)
\(\gamma\)	Scaling exponent for the effect of \(\lambda\) on \(\Delta t_r\)

---

Manuscript Figures

# rescaled graphs
legend <- get_legend(heating_graph + theme(legend.box.margin = margin(0, 0, 0, 0)))

p2 <- plot_grid(heating_graph + theme(legend.position="none"),
                rescaled_duration + theme(legend.position="none"),
                # rescaled_tolerances + theme(legend.position="none"),
                labels = "AUTO", nrow = 1, align = "hv", hjust = -1, label_size = 10)

fig2 <- plot_grid(p2, legend, rel_widths =c(3, 0.5))
fig2

# prediction graphs
pred_high_low <- pred3 + aes(col = heating_rate) + 
  scale_x_continuous(expand = c(0,0), limits = c(0, 40)) +
  scale_y_continuous(expand = c(0,0), limits = c(0, 40)) +
  geom_smooth(method = "lm", se = FALSE, col = 1) +
  scale_colour_viridis_c(name = expression("Heating rate (\u00B0C h" ^-1 * ")"), trans = "log10", labels = function(x) signif(x)) +
  ggtitle("Predict from lower acclimation temperature") +
  theme(legend.position="none")

pred_low_high <- pred4 + aes(col = heating_rate) + 
  scale_x_continuous(expand = c(0,0), limits = c(0, 40)) +
  scale_y_continuous(expand = c(0,0), limits = c(0, 40)) +
  geom_smooth(method = "lm", se = FALSE, col = 1) +
  scale_colour_viridis_c(name = expression("Heating rate (\u00B0C h" ^-1 * ")"), trans = "log10", labels = function(x) signif(x)) +
  ggtitle("Predict from higher acclimation temperature") +
  theme(legend.position="none")

p3 <- plot_grid(pred_high_low, pred_low_high, labels = "AUTO", nrow = 1, align = "h", hjust = -1, label_size = 10)
fig3 <- plot_grid(p3, legend, rel_widths =c(3, 0.5))
fig3

# Heating rate effect
fig5 <- plot_grid(g_departure, g_scaling, g_predict_tau, labels = "AUTO", nrow = 1, align = "h", hjust = -1, label_size = 10)
fig5

Figures saved as TIFF. Resulting dataset (ARR) saved as CSV.

---

Session info

R version 4.2.2 (2022-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19044)

Matrix products: default

locale:
[1] LC_COLLATE=English_United Kingdom.utf8 
[2] LC_CTYPE=English_United Kingdom.utf8   
[3] LC_MONETARY=English_United Kingdom.utf8
[4] LC_NUMERIC=C                           
[5] LC_TIME=English_United Kingdom.utf8    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] cowplot_1.1.1   forcats_0.5.2   stringr_1.4.1   dplyr_1.0.10   
 [5] purrr_0.3.5     readr_2.1.3     tidyr_1.2.1     tibble_3.1.8   
 [9] ggplot2_3.4.0   tidyverse_1.3.2

loaded via a namespace (and not attached):
 [1] lattice_0.20-45     lubridate_1.9.0     assertthat_0.2.1   
 [4] rprojroot_2.0.3     digest_0.6.30       utf8_1.2.2         
 [7] R6_2.5.1            cellranger_1.1.0    backports_1.4.1    
[10] reprex_2.0.2        evaluate_0.18       highr_0.9          
[13] httr_1.4.4          pillar_1.8.1        rlang_1.0.6        
[16] googlesheets4_1.0.1 readxl_1.4.1        rstudioapi_0.14    
[19] jquerylib_0.1.4     Matrix_1.5-1        rmarkdown_2.18     
[22] labeling_0.4.2      splines_4.2.2       googledrive_2.0.0  
[25] munsell_0.5.0       broom_1.0.1         compiler_4.2.2     
[28] modelr_0.1.9        xfun_0.34           pkgconfig_2.0.3    
[31] mgcv_1.8-41         htmltools_0.5.3     tidyselect_1.2.0   
[34] bookdown_0.30       viridisLite_0.4.1   fansi_1.0.3        
[37] crayon_1.5.2        tzdb_0.3.0          dbplyr_2.2.1       
[40] withr_2.5.0         grid_4.2.2          nlme_3.1-160       
[43] jsonlite_1.8.3      gtable_0.3.1        lifecycle_1.0.3    
[46] DBI_1.1.3           magrittr_2.0.3      scales_1.2.1       
[49] cli_3.4.1           stringi_1.7.8       cachem_1.0.6       
[52] farver_2.1.1        fs_1.5.2            xml2_1.3.3         
[55] bslib_0.4.1         ellipsis_0.3.2      generics_0.1.3     
[58] vctrs_0.5.0         tools_4.2.2         glue_1.6.2         
[61] hms_1.1.2           fastmap_1.1.0       yaml_2.3.6         
[64] timechange_0.1.1    colorspace_2.0-3    gargle_1.2.1       
[67] rvest_1.0.3         knitr_1.40          haven_2.5.1        
[70] sass_0.4.2         

---

References

Dell, Anthony I., Samraat Pawar, and Van M. Savage. 2011. “Systematic Variation in the Temperature Dependence of Physiological and Ecological Traits.” Journal Article. Proceedings of the National Academy of Sciences 108 (26): 10591–96. https://doi.org/10.1073/pnas.1015178108.

Morley, S. A., A. E. Bates, M. Lamare, J. Richard, K. D. Nguyen, J. Brown, and L. S. Peck. 2016. “Rates of Warming and the Global Sensitivity of Shallow Water Marine Invertebrates to Elevated Temperature.” Journal Article. Journal of the Marine Biological Association of the United Kingdom 96 (1): 159–65. https://doi.org/10.1017/S0025315414000307.

Morley, S. A., L. S. Peck, J. M. Sunday, S. Heiser, and A. E. Bates. 2019. “Physiological Acclimation and Persistence of Ectothermic Species Under Extreme Heat Events.” Journal Article. Global Ecology and Biogeography 28 (7): 1018–37. https://doi.org/10.1111/geb.12911.

Morley, S. A., L. S. Peck, J. Sunday, S. Heiser, and A. E. Bates. 2018. “Acclimation Potential of Global Ectothermic Species, Collated from Literature, 1960 to 2015.” Dataset. Natural Environment Research Council, Cambridge, UK. https://doi.org/https://doi.org/10.5285/20010bfb-c6d3-430f-b1f7-d16790ab8359.

R Core Team. 2022. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.

Wickham, Hadley. 2019. Tidyverse: Easily Install and Load the Tidyverse. https://CRAN.R-project.org/package=tidyverse.

Wilke, Claus O. 2020. Cowplot: Streamlined Plot Theme and Plot Annotations for ’Ggplot2’. https://CRAN.R-project.org/package=cowplot.

Xie, Yihui. 2020. Bookdown: Authoring Books and Technical Documents with r Markdown. https://github.com/rstudio/bookdown.

---

1. Trinity College Dublin, Ireland, kongj@tcd.ie↩︎

2. Trinity College Dublin, Ireland, jacksoan@tcd.ie↩︎

3. Trinity College Dublin, Ireland, paynen@tcd.ie↩︎

4. CNRS, France, arnoldi.jeff@gmail.com↩︎