\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[ngerman,greek,english]{babel}



\usepackage{float}








\begin{document}

\title{Variation in the rate of phenotypic plasticity among ectotherms}



\author[1]{Sigurd Einum}%
\author[1]{Tim Burton}%
\affil[1]{Norwegian University of Science and Technology}%


\vspace{-1em}



  
  \date{June 16, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
An individual's fitness cost associated with environmental change
depends on the rate of adaptive phenotypic plasticity, and yet our
understanding of plasticity rates in an ecological and evolutionary
context remains limited. We provide the first quantitative synthesis of
existing plasticity rate data, focusing on acclimation of temperature
tolerance in ectothermic animals, where we demonstrate applicability of
a recently proposed analytical approach. The analyses reveal
considerable variation in plasticity rates of this trait among species,
with half-times (how long it takes for the initial deviation from the
acclimated phenotype to be reduced by 50\% when individuals are shifted
to a new environment) ranging from 3.7 to 770.2 h. Furthermore, rates
differ among higher taxa, being higher for amphibians and reptiles than
for crustaceans and fishes, and with insects being intermediate. We
argue that a more comprehensive understanding of phenotypic plasticity
will be attained through increased focus on the rate parameter.%
\end{abstract}%



\sloppy


\textbf{Variation in the rate of phenotypic plasticity among ectotherms}

\par\null

Sigurd Einum\textsuperscript{1*}and Tim Burton\textsuperscript{1,2}

\textsuperscript{1}Centre for Biodiversity Dynamics, Department of
Biology, Norwegian University of Science and Technology, Realfagbygget,
NO-7491 Trondheim, Norway

\textsuperscript{2}Norwegian Institute for Nature Research, PO Box 5685,
Torgarden, NO-7485, Trondheim, Norway.

*corresponding author: sigurd.einum@ntnu.no

Type of article: Letter

Number of words in abstract: 149

Number of words in main text: 4500

Number of references: 30

Number of figures: 2

Number of tables: 2

Statement of authorship: both authors conceived the study and identified
data sources, SE extracted the data, conducted the analyses, and wrote
the first draft of the manuscript, both authors contributed to
revisions.

Data accessibility statement: The data will be published in Dryad upon
acceptance

Running title: Variation in phenotypic plasticity rates

Keywords: Phenotypic plasticity, acclimation rate, comparative,
temperature, time course, body size, selection, fitness, exponential
decay, critical



\textbf{Abstract}

An individual's fitness cost associated with environmental change
depends on the rate of adaptive phenotypic plasticity, and yet our
understanding of plasticity rates in an ecological and evolutionary
context remains limited. We provide the first quantitative synthesis of
existing plasticity rate data, focusing on acclimation of temperature
tolerance in ectothermic animals, where we demonstrate applicability of
a recently proposed analytical approach. The analyses reveal
considerable variation in plasticity rates of this trait among species,
with half-times (how long it takes for the initial deviation from the
acclimated phenotype to be reduced by 50\% when individuals are shifted
to a new environment) ranging from 3.7 to 770.2 h. Furthermore, rates
differ among higher taxa, being higher for amphibians and reptiles than
for crustaceans and fishes, and with insects being intermediate. We
argue that a more comprehensive understanding of phenotypic plasticity
will be attained through increased focus on the rate parameter.



\textbf{Introduction}

Phenotypic plasticity, or the expression of different phenotypes across
environments by a single genotype, is an important process by which
organisms can minimize environmental impacts on fitness (Gabriel, 2005;
Gabriel, Luttbeg, Sih, \& Tollrian, 2005; Padilla \& Adolph, 1996;
Siljestam \& \selectlanguage{ngerman}Östman, 2017). Such plasticity can be described by two
parameters. First, the capacity for plasticity determines the amount by
which the phenotype can change following a shift in the environment.
This parameter can be measured as the change in the slope of the
relationship between trait value and environment as plasticity proceeds,
from acute exposure until the full plastic response has been achieved
(see Einum et al. 2019 for arguments why it is this change in the slope,
and not slope \emph{per se} , that describes plasticity). Whereas the
capacity for plasticity has received considerable theoretical and
empirical interest from ecologists and evolutionary biologists,
empirical support for certain predictions regarding the evolution of
this plasticity parameter remain equivocal. For example, it has been
proposed that organisms inhabiting more variable environments should
evolve greater plasticity capacities. However, this is rarely supported
by empirical data (Gunderson \& Stillman, 2015; Kelly, Sanford, \&
Grosberg, 2012; MacLean et al., 2019; Pereira, Sasaki, \& Burton, 2017;
Phillips et al., 2016; Sgro et al., 2010; van Heerwaarden, Kellermann,
\& Sgrò, 2016; van Heerwaarden, Lee, Overgaard, \& Sgrò, 2014).
Recently, Burton et al. (2022) suggested that this discrepancy gives
reason for pause, and that greater considerations of the second
parameter, the rate of plasticity, which addresses the timescale over
which plastic phenotypic change occurs, might aid in bringing this field
of research forward.

If the plasticity of a trait is an adaptive response, the fitness cost
that an organism incurs following a change in its environment should be
minimized once the phenotype becomes fully adjusted to the new
environment. Hence, the rate at which the phenotype approaches this
state should determine how long the individual expresses a sub-optimal
phenotype, and in part, determine the magnitude of the fitness cost
associated with that change in the environment. Given that organisms are
unlikely to be able to predict changes in all of the relevant
environmental variables they are exposed to, it seems plausible that
individuals may actually spend a considerable proportion of their time
having a phenotype that is not fully adjusted to their current
environment. This mismatch between environment and phenotype, and
associated cumulative fitness costs, will be exacerbated if plastic
responses are slow relative to the timescale of environmental change.
Furthermore, as pointed out by Burton et al. (2022), the rate of
plasticity might even influence how the capacity for plasticity evolves,
because the evolution of capacity depends on the predictability of the
environment. Organisms that can rapidly implement their phenotypic
response to a new environment can postpone the onset of this process
closer to the time of selection in that environment than organisms that
do so at a slower rate. In a temporally autocorrelated environment this
would in effect make faster responding organism to more accurately
`predict' the future selective environment at the moment when they have
to start adjust their phenotype. Thus, a faster rate of plasticity might
effectively increase predictability in the environment, which in turn
should favour the evolution of greater phenotypic plasticity capacity
(Lande 2014).

Presently, a basic synthesis of plasticity rate data among taxa is
lacking, and consideration of how rates of plasticity might be expected
to evolve in response to environmental change is absent from current
theoretical models (Lande 2014, Siljestam \& Östman 2017). This
knowledge gap is mirrored by, and perhaps stems from, a lack of
empirical interest in studying the evolution of rates of plasticity.
Although a substantial number of empirical studies document how
phenotypes change over time when introduced into new environments, these
studies remain largely descriptive, fail to address evolutionary
hypotheses, and very rarely (four out of 166 studies surveyed by Burton
et al. 2022) attempt to provide any formal statistical quantification of
the time course of plasticity. Thus, advancing our understanding of the
evolution of phenotypic plasticity might arguably benefit from a shift
in focus from capacities to rates of plasticity. To stimulate such a
shift, we provide the first comparative analysis of published data
describing rates of plasticity. In doing so, we follow recent
suggestions (Burton et al. 2022) regarding the estimation of plasticity
rates in a (i) standardized way, which is (ii) consistent with theory
and (iii) directly comparable across taxa and traits.

We draw upon published data from studies of acclimation to temperature
among ectotherms. Temperature is an environmental variable that affects
all organisms, varies substantially in space and time, and which has
particularly pervasive effects on biochemical, physiological and
ecological processes in this group of animals (Daufresne et al. 2009).
We focus our synthesis on traits describing temperature tolerance. We
first determine the shape of how temperature tolerance changes over time
(exponential vs. linear decay) in response to a shift in ambient
temperature, as this is the first step required when calculating the
rate of plasticity. After calculating rates of plasticity for each
published dataset, we then investigate relationships between rates of
plasticity and taxonomic class, body size, and acclimation temperature.
By providing clear evidence that rates of plasticity have diverged among
ectotherm classes we show that the rate of plasticity can, and does,
evolve, and that increased empirical and theoretical focus on the rate
parameter is likely to provide a way forward in understanding evolution
of phenotypic plasticity.

\textbf{Methods}

\emph{Collection of data}

We identified published studies documenting the rate of change in
whole-body measurements of thermal tolerance in response to manipulation
of ambient temperature. Specifically, we sought to identify papers where
the study organisms, having been acclimated to an initial temperature
for a given period of time, were transferred to a `new' temperature,
after which thermal tolerance was measured on subsamples of these
individuals at different points in time. Our search criteria included
studies that addressed tolerance of both high and low temperatures and
performed measurements of critical temperatures (i.e. CTmax, CTmin),
time to death or immobility, or mortality rate at a stressful
temperature. An additional criterion was that measurements of thermal
tolerance were made at a minimum of three time points. Using these
criteria, relevant papers were identified using the following procedure.
First, we included eight papers identified by Burton et al. (2022) that
were based on a search in Web of Science (see Burton et al. for
procedure). Second, we conducted two searches on Google Scholar during
January 2022. In the first search, we used the terms ''temperature'' and
''rate of acclimation''. The 739 hits were sorted based on relevance,
and the first 200 abstracts were scanned for inclusion, resulting in a
total of 55 papers (none of these appeared among the last 30 scanned
abstracts). In the second search, we used the terms ''temperature'',
''acclimation'', ''time course'' and ''critical''. This resulted in
30200 hits, which also were sorted based on relevance, and where the
first 200 abstracts were scanned for inclusion. For this second search,
7 additional papers were identified (the majority of relevant hits had
been identified in the initial search). Finally, the reference lists of
papers identified above were scanned for reference to additional papers
not located in any of the previous searches, resulting in a further 18
papers. Thus, a total of 88 papers were identified during this
procedure. For each paper, we determined whether it presented the
required data in an accessible form. Fifty-nine of the papers identified
above were found to contain suitable data. Many of these contained data
from several experiments (different species and/or different acclimation
temperatures), and data from a total of 308 acclimation experiments from
a total of 114 species were extracted from these papers. For most of the
papers, data were presented in figure format only, usually as the mean
tolerance of groups of individuals that had experienced the `new'
temperature for differing lengths of time prior to being measured for
the chosen tolerance trait. We digitized such data using the
WebPlotDigitizer (https://apps.automeris.io/wpd/). In addition to
thermal tolerance data we extracted information for the following
variables: taxonomic class, species name, mean body mass of individuals
used in the experiment, and acclimation temperature. For experiments
where body mass information was absent we used the following procedure
to obtain such data. First, we searched for information on minimum and
maximum adult body size (length or mass) for the species in question. If
only data on body length were available, we searched for allometries
between length and mass for that species, or for species within the same
genus, and used these to convert body lengths to mass. We then
calculated the mean of the minimum and maximum body mass values. If no
data on minimum or maximum size was found, we searched for studies that
presented body sizes of adult individuals during field surveys and used
the mean of these. Combining body size estimates from each of these
sources yielded measurements for 290 of the 308 experiments described
above, representing a total of 107 species from the following classes of
ectotherms; amphibians, reptiles, insects, osteichtyes and malacostraca.
For simplicity the latter two will be referred to as fishes and
crustaceans, respectively.

\emph{Calculating rates of plasticity}

Rates of plasticity were calculated according to the approach outlined
by Burton et al. (2022). Thus, for each observation of temperature
tolerance at times \emph{t} we first calculated the proportion of the
full plastic response that remained to be
achieved,\emph{D\textsubscript{t}} , as (\emph{z\textsubscript{t}}
-\selectlanguage{english}\emph{z\textsubscript{[?]}} )/(\emph{z\textsubscript{0}}
--\selectlanguage{english}\emph{z\textsubscript{[?]}} ), where\emph{z\textsubscript{0}} is the
first measurement of the phenotype (typically measured prior to the
onset of acclimation), \emph{z\textsubscript{t}} is the phenotype after
acclimating for a time period \emph{t} , and \selectlanguage{english}\emph{z\textsubscript{[?]}}
is the fully adjusted phenotype. In theory, \selectlanguage{english}\emph{z\textsubscript{[?]}} is
the final measurement of thermal tolerance. However, given measurement
noise this may not be the case. Thus, we
defined\selectlanguage{english}\emph{z\textsubscript{[?]}} as the maximum (for acclimation to
higher temperature) or minimum (for acclimation to lower temperature)
value for a given measurement of tolerance observed. For each
experiment, we considered two potential responses in plasticity of
thermal tolerance. First, \emph{D\textsubscript{t}} may have a value of
1 at \emph{t} = 0 (first measurement) and decline linearly with time
with a rate\selectlanguage{greek}\emph{λ\textsubscript{L}}\selectlanguage{english} given by \selectlanguage{greek}\emph{λ\textsubscript{L}}\selectlanguage{english}
= (1-\emph{D\textsubscript{t}} )/\emph{t} until it reaches 0, after
which it will be constant. Alternatively, \emph{D\textsubscript{t}} may
have a value of 1 at \emph{t} = 0 and decline as an exponential decay
function towards 0 with rate \selectlanguage{greek}\emph{λ\textsubscript{E}}\selectlanguage{english} , such
that\(D_{t}=e^{{-\lambda}_{E}t}\). Thus, for each experiment, we fitted these two
types of models to the observed data. For the linear decline we fitted a
piecewise regression with an intercept of 1 (at \emph{t} = 0), which
estimates the breakpoint\emph{b} (i.e. the time at which
\emph{D\textsubscript{t}} reaches 0), from which the rate
\selectlanguage{greek}\emph{λ\textsubscript{L}}\selectlanguage{english} can be calculated as 1/\emph{b} . Thus, both
models estimate a single parameter, and their relative fits for a given
experiment can be compared directly using their residual standard
errors. Models were fitted using the function\emph{nls\_multstart} from
the \emph{nls.multstart} package (Padfield and Matheson 2020).

For most studies, \emph{z\textsubscript{0}} was measured in individuals
prior to transfer to the new temperature. It was often less clear
whether experimenters had been able to measure a `true' value
of\emph{z\textsubscript{[?]}} , i.e. thermal tolerance after full
acclimation to the new temperature had been obtained. This is a key
point when measuring rates of plasticity because estimates
of\selectlanguage{greek}\emph{λ\textsubscript{E}}\selectlanguage{english} will be biased if full acclimation to the
new environment has not been achieved (Fig. S1). However, an advantage
of the exponential decay function is that achievement of full
acclimation can be assessed by calculating the slope of the estimated
function (i.e.\({{-\lambda}_{E}e}^{\lambda_{E}t}\)) at the final acclimation time point
\emph{t\textsubscript{n}} (this slope has an asymptotic value of 0).
Thus, this value was included as a covariate in our analysis (see below)
to control for any bias introduced by variation in maximum acclimation
time among studies.

\emph{Statistical analyses}

All statistical analyses were conducted in R v.4.1.2 (R Core Team,
2021). Comparisons of exponential and linear decay models showed that
the former was superior to the latter in most cases (see Results). Thus,
the following statistical analyses use \selectlanguage{greek}\emph{λ\textsubscript{E}}\selectlanguage{english} as the
dependent variable.

Inspection of the distribution of standard errors (SE) for estimates
of\selectlanguage{greek}\emph{λ\textsubscript{E}}\selectlanguage{english} were used to exclude experiments that had
poor model fit for the exponential decline models (Fig. S2). This was
based on judgement, trading-off maintaining sample size while avoiding
inclusion of potentially biased estimates. Using a threshold SE value of
0.01 enabled us to retain 78\% of the experiments (n = 240 out of 307).
The statistical analyses described below use estimates obtained below
this threshold. Increasing the threshold SE value to 0.02 increased
inclusion rate to 88\% (n = 271) without resulting in qualitative
changes in the results (Table S1). The bias in
estimated\selectlanguage{greek}\emph{λ\textsubscript{E}}\selectlanguage{english} for experiments with SE
\textgreater{} 0.02 and within the range 0.01-0.02 can be observed by
comparing\selectlanguage{greek}\emph{λ\textsubscript{E}}\selectlanguage{english} in experiments within different SE
bins (Fig. S3).

Variation in \selectlanguage{greek}\emph{λ\textsubscript{E}}\selectlanguage{english} among taxonomic classes of
ectothermic animals was analysed using linear mixed-effect models,
fitted using the function \emph{lme} in the package \emph{nlme}(Pinheiro
et al., 2022). The full model contained the fixed effects of class, body
mass, acclimation temperature, and slope of the estimated exponential
decay function at \emph{t\textsubscript{n}} , whereas species was
included as a random effect.

The full model was compared to simplified ones based on AICc values
using the function \emph{dredge} from the package \emph{MuMIn} (Barton
2020). For all analyses inspection of residual plots suggested that
assumptions of their normality and homogeneity were satisfied.

\textbf{Results}

A total of 290 experiments could be fitted with both piecewise
regression and exponential decay models (Fig. S4). The remaining 18
could not be fitted with piecewise regression, which most frequently
occurred because the minimum value of \emph{D} (i.e. complete
acclimation) was reached at the second observation, in which case
multiple breakpoints can yield an identical fit. Of these 290
experiments, comparisons of residual standard errors from models showed
that the data in 62\% (n=180) of them were best explained by exponential
decay models. Furthermore, superior fit by piecewise regressions (N =
110 experiments) was observed primarily in experiments where thermal
tolerance was measured at relatively few time points following exposure
of the study organisms to the new temperature (Fig. S5). Thus, for more
comprehensive experiments (i.e. greater number of measurement points in
time), the shape of the plasticity rate response is better described by
an exponential decay function. Thus, in all further analyses we use the
estimates of \selectlanguage{greek}\emph{λ\textsubscript{E}}\selectlanguage{english} (Table S2).

For the data included in the statistical analyses presented here,
species-specific mean estimated \selectlanguage{greek}\emph{λ\textsubscript{E}}\selectlanguage{english} was 0.0342
h\textsuperscript{-1}. The variation in estimated
mean\selectlanguage{greek}\emph{λ\textsubscript{E}}\selectlanguage{english} across species was considerable, with SD =
0.0305 (i.e. CV = 89\%). Species-specific mean
\selectlanguage{greek}\emph{λ\textsubscript{E}}\selectlanguage{english}ranged from 0.0009 to 0.1892
h\textsuperscript{-1}. Variation in\selectlanguage{greek}\emph{λ\textsubscript{E}}\selectlanguage{english} was best
explained by a model that included taxonomic class, acclimation
temperature, and the slope of the estimated exponential decay function
at \emph{t\textsubscript{n}} , while evidence for an effect of body mass
was weak (Table 1). Model coefficients of the best ranked model show
that plasticity rates were highest among amphibians and reptiles, lowest
among fishes and crustaceans, and intermediate in insects (Table 2).
This pattern was also evident when examining distributions of
\selectlanguage{greek}\emph{λ\textsubscript{E}}\selectlanguage{english} among these classes without correcting for the
covariates fitted in the best ranked model (Fig. 1). We also observed
that plasticity rates increased with acclimation temperature (Table 1,
2, Fig. 2). Finally, plasticity rates were observed to be higher when
the slope of the estimated exponential decay function at
\emph{t\textsubscript{n}} was shallower (Table 1). In other words,
experiments where complete acclimation was more likely to have been
obtained were associated with higher estimated plasticity rates. This
pattern was mainly driven by the amphibian data, which had a large
number of experiments with relatively steep slopes
at\emph{t\textsubscript{n}} (Fig. S6). This is in the opposite direction
to what might be predicted if a relationship between this slope and
estimated acclimation rate is a statistical artefact (Fig. S1). Rather,
the causality of this relationship is likely in the opposite direction,
i.e. experiments that had a low acclimation rate were
stopped\emph{before} complete acclimation to the new temperature had
occurred. Nevertheless, to evaluate if this anomaly could have
influenced our results, we repeated the model comparison summarized in
Table 2 but including only experiments where the slope of the estimated
exponential decay function at \emph{t\textsubscript{n}} was larger than
-0.002 and -0.001, respectively. For both cases, results remained
qualitatively identical and showed only minor quantitative changes
(Table S3 - S6).

\textbf{Discussion}

Here, we provide the first synthesis of plasticity rates (\selectlanguage{greek}\emph{λ} \selectlanguage{english})
across taxa and apply a novel method to obtain standardized and thus
comparable measures of this plasticity parameter. For our focal trait,
acclimation of temperature tolerance to temperature change, the shape of
the plasticity responses were well described by an exponential decay
function (i.e. with the rate \selectlanguage{greek}\emph{λ\textsubscript{E}}\selectlanguage{english} ). In other
words, the absolute rate of change in temperature tolerance when an
individual is shifted to a new temperature is proportional to the
deviance from the phenotype when completely acclimated to that
temperature. We thus validate an assumption that has previously been
made in theory describing the evolution of phenotypic plasticity (Lande
2014). In contrast, superior fit by piecewise linear regression was
primarily observed in experiments with poor temporal resolution, thus
demonstrating the importance of measuring phenotypes at multiple time
points in the new environment to determine the optimal model to use when
estimating the shape of the plasticity response.

Variation in estimated \selectlanguage{greek}\emph{λ\textsubscript{E}}\selectlanguage{english} among species was
considerable. To put these numbers into perspective they can be
translated into half-times, or how long it takes for the initial
deviation from the fully acclimated phenotype to be reduced by 50\%
after being shifted to a new environment, which is given as
ln(2)/\selectlanguage{greek}\emph{λ\textsubscript{E}}\selectlanguage{english} . The mean
observed\selectlanguage{greek}\emph{λ\textsubscript{E}}\selectlanguage{english} of 0.0342 h\textsuperscript{-1}
corresponds to a half-time of 20.3 h, whereas the minimum and maximum
species-specific\selectlanguage{greek}\emph{λ\textsubscript{E}}\selectlanguage{english} , estimated to 0.0009 and
0.1892 h\textsuperscript{-1}, correspond to half-times of 770.2 and 3.7
h, respectively (only including data used in the analyses, i.e. with
SE\selectlanguage{greek}\emph{λ\textsubscript{E}}\selectlanguage{english} \textless{} 0.01). Using model estimates
(Table 2), half-times at 20 \selectlanguage{ngerman}°C when the final slope of the decay
function is zero is 10.7, 15.8, 30.1, 43.3 and 40.8 h for amphibians,
reptiles, insects, crustaceans and fishes, respectively. Thus, our
analyses demonstrate considerable systematic variation in rates of
phenotypic plasticity among these taxa. This begs the question of why
such variation has evolved. One might speculate that due to the higher
heat capacity of water compared to air, fishes (all species) and
crustaceans (5 out of 7 species in our data) which inhabit aquatic
environments experience higher temperature stability than taxa
inhabiting terrestrial environments. Furthermore, aquatic habitats are
often thermally stratified (lentic habitats such as lakes and oceans) or
provide cold-water plumes (rivers), which allows for behavioural
thermoregulation under periods of stressful temperatures (Kurylyk et al.
2015, Freitas et al. 2016, Harrison et al. 2016). This may reduce the
strength of selection on rapid plasticity compared to in the terrestrial
environment occupied by reptiles and insects (all species in our data).
Although the amphibians also inhabit aquatic environments (particularly
during juvenile life stages), their utilization of thermal refugia in
deep or fast flowing water is likely limited. Unfortunately, habitat use
is confounded with phylogeny in this data set, preventing direct
analysis of the effect of habitat on rate of plasticity in thermal
tolerance. Thus, although we do demonstrate evolutionary divergence in
plasticity rates among these taxa, it remains an open question as to
whether this pattern results from evolutionary adaptation to
environmental conditions. However, this question could be addressed in
future work that targets populations or species that experience known
and contrasting patterns of environmental variability.

We observed a positive relationship between acclimation temperature and
rate of plasticity in temperature tolerance. This pattern may be
explained by the general relationship that exists between developmental
rate and body temperature in ectotherms, which is driven by the positive
effect of temperature on biochemical reactions and metabolic rate (Brown
et al. 2004). It may also explain the observation that within a species,
acclimation to high temperature is achieved faster than acclimation to
low temperature (Burton et al. 2020). A relationship between metabolic
rate and the rate of plasticity was previously hypothesized and
addressed by Rohr et al. (2018), but in a less direct manner.
Specifically, Rohr et al. (2018) argued that the effect of metabolic
rate on rates of thermal plasticity in ectotherms should be evident as a
negative relationship between body size and plasticity rate, because
smaller organisms tend to have a higher mass-specific metabolic rate
than larger ones. They did not however, calculate rates of plasticity
from experiments that were explicitly designed to do so. Rather, they
used data from experiments that measure the phenotype at only two time
points (\emph{z\textsubscript{0}} and \selectlanguage{english}\emph{z\textsubscript{[?]}} in our
terminology), and from this inferred how the bias in acclimation
capacity caused by insufficient acclimation time was influenced by body
size. Based on their results it was concluded that rates of plasticity
appeared to be higher for smaller organisms. Using a more direct
approach we failed to find support for a general relationship between
body size and rate of thermal plasticity. Yet, our observation that the
rate of plasticity in temperature tolerance is positively related to
acclimation temperature suggests a role for metabolic rate in causing
some of the variation in plasticity rate across experiments.

Given the general patterns in rate of plasticity observed here, further
efforts in studying this plasticity parameter may be fruitful and
provide a better foundation for understanding how plasticity evolves in
response to environmental variation. From an empirical perspective,
including a temporal-dimension in experiments that study plasticity may
be included without large costs. In this respect, we make two
recommendations. First, a proper choice of model (linear vs. exponential
decay) for estimating lambda requires multiple measurements of the
phenotype as it responds to the new environment. Our analyses indicate
that five or more measurements may be required to adequately establish
the shape of the plasticity response (Fig. S5). Superficially, this
requirement might appear to substantially increase the workload of such
studies in comparison to studies that only estimate the capacity for
plastic phenotypic change in a trait. However, once the model that best
describes the shape of the plastic response to the new environment is
established, a single measurement \emph{z\textsubscript{t}} after
time\emph{t} (which must be prior to achievement of full acclimation) in
addition to those typically measured (\emph{z\textsubscript{0}} in
non-acclimated individuals and \selectlanguage{english}\emph{z\textsubscript{[?]}} after the full
acclimation response has been obtained) is sufficient to accurately
estimate \emph{D\textsubscript{t}} , which in turn can be used to
calculate the rate of plasticity (\selectlanguage{greek}\emph{λ\textsubscript{E}}\selectlanguage{english} =
ln(\emph{D\textsubscript{t}} )/\emph{t} for exponential decay
or\selectlanguage{greek}\emph{λ\textsubscript{L}}\selectlanguage{english} = (1-\emph{D\textsubscript{t}} )/\emph{t}
for linear decay). Thus, the workload in such experiments can be greatly
reduced by performing a pilot experiment with sufficient temporal
resolution (in terms of measurement time points) that provides a precise
description of the shape of plasticity response to the new environment
before performing more replicates at a lower temporal resolution to
obtain the desired estimates of \selectlanguage{greek}\emph{λ} . \selectlanguage{english}It should be noted
that\selectlanguage{greek}\emph{λ\textsubscript{E}}\selectlanguage{english} and \selectlanguage{greek}\emph{λ\textsubscript{L}}\selectlanguage{english} are not
directly comparable, because the initial approach towards the fully
adjusted phenotype is more rapid under exponential decay. Thus, the
relative support for these two types of plasticity responses should be
reported. As a second recommendation, experimenters should strive to
ensure that complete acclimation to the new environment is achieved
prior to measuring \emph{z\textsubscript{[?]}} . Our analyses show that
failing to do so can, and does, lead to bias in estimation
of\selectlanguage{greek}\emph{λ\textsubscript{E}}\selectlanguage{english} (Fig. S1, S6). Ideally this is achieved by
rearing individuals in all the alternative environments for the whole
duration of the experiment (i.e. both prior to and after some of the
individuals are transferred into new environments). As pointed out by
Burton et al. (2022), this has rarely been done in studies of rates of
plasticity. Rather, the majority of studies first acclimate the animals
to a single initial environment before shifting them to a new
environment and then performing repeated measures of phenotype in this
new environment for what typically appears to be a pre-determined (and
potentially insufficient) duration.

Natural next steps in research on evolution of plasticity would be to
test for links between environmental variation and the evolution of
rates of plasticity, and to provide theoretical models that address the
co-evolution of plasticity rates and capacity (see Introduction).
Although this is beyond the scope of the current paper, our work
provides both methodology and novel insights that should stimulate
future work along these lines. We also re-emphasize a point made
previously (Burton et al. 2022) - that selection on the rate of
plasticity might be stronger than selection on the capacity for
plasticity. Evolutionary theory posits a central role for phenotypic
plasticity in mitigating the fitness impact of environmental variation,
but that possessing the potential for such a response is associated with
a fitness cost in stable environments (Lande 2009). Fitness costs of
plasticity can be categorized into costs of maintenance and costs of
production. Costs of maintenance represent the investment of resources
into maintaining the machinery required for detecting and responding to
a change in the environment and will be paid at a constant rate
independent of environmental conditions (Auld et al. 2010). In contrast,
production costs are only paid when the plastic response is triggered
and are compensated by the fitness benefits associated with changing the
phenotype. If one assumes that the capacity for plasticity can be
increased by operating the `machinery' required to change a trait for a
longer duration, this will increase production costs but not maintenance
costs. Populations living in less variable environments may therefore
pay a small price for maintaining their capacity for plasticity (as
shown by Van Buskirk \& Steiner 2009), and adaptation of this parameter
of plasticity to levels of environmental fluctuations may therefore be
relatively modest in magnitude. In contrast, increasing the rate of
change in the same trait would require increasing the size or output of
that `machinery', with corresponding increases in maintenance costs.
Populations living in less variable environments should therefore
experience strong selection against maintaining rapid plasticity due to
higher maintenance costs, and adaptative evolution across populations
may then be expected to be more pronounced for the rate of plasticity.
This line of reasoning is also consistent with theoretical results
showing that maintenance costs shape the evolution of plasticity to a
greater extent than production costs (Sultan \& Spencer 2002). Given
these considerations, and the results presented in the current study, it
seems prudent to address the hypothesis that adaptation to environmental
variation may be more pronounced in terms of rates of plasticity rather
than capacity of plasticity. By providing clear evidence that rates of
plasticity have diverged among ectotherm classes we show that it is a
trait that evolves, and that increased empirical and theoretical focus
on the rate parameter is likely to provide a way forward for a more
comprehensive understanding of phenotypic plasticity.

\textbf{References}

Auld, J. R., Agrawal, A. A., \& Relyea, R. A. (2010). Re-evaluating the
costs and limits of adaptive phenotypic plasticity. \emph{Proceedings of
the Royal Society B: Biological Sciences, 277} (1681), 503-511.

Barton, K. (2020). MuMIn: Multi-model inference. R package version
1.43.17. https://CRAN.R-project.org/package=MuMIn

Brown, J.H., Gillooly, J.F., Allen, A.P., Savage, V.M. and West, G.B.
(2004). Toward a metabolic theory of ecology. Ecology, 85:
1771-1789.~https://doi.org/10.1890/03-9000

Burton T., Lakka H.-K., Einum S. (2020). Acclimation capacity and rate
change through life in the zooplankton \emph{Daphnia} . Proc. R. Soc. B
287: 20200189. http://dx.doi.org/10.1098/rspb.2020.0189

Burton, T., Ratikainen, I. I., Einum, S. (2022). Environmental change
and the rate of phenotypic plasticity. Glob Change Biol, in press.
Preprint: https://doi.org/10.6084/m9.figshare.18741044.v1

Daufresne, M., Lengfellner, K., \& Sommer, U. (2009). Global warming
benefits the small in aquatic ecosystems. \emph{Proc Natl Acad Sci, 106}
, 12788-12793. https://doi.org/10.1073/pnas.0902080106

Einum, S., Ratikainen, I., Wright, J., P\selectlanguage{ngerman}élabon, C., Bech, C., Jutfelt,
F., Stawski, C. and Burton, T. (2019), How to quantify thermal
acclimation capacity?. Glob Change Biol, 25:
1893-1894.~https://doi.org/10.1111/gcb.14598

Freitas, C., Olsen, E.M., Knutsen, H., Albretsen, J. and Moland, E.
(2016), Temperature-associated habitat selection in a cold-water marine
fish. J Anim Ecol, 85: 628-637.~https://doi.org/10.1111/1365-2656.12458

Gabriel, W. (2005). How stress selects for reversible phenotypic
plasticity. \emph{Journal of Evolutionary Biology, 18} (4), 873-883.
doi:10.1111/j.1420-9101.2005.00959.x

Gabriel, W., Luttbeg, B., Sih, A., \& Tollrian, R. (2005). Environmental
tolerance, heterogeneity, and the evolution of reversible plastic
responses. \emph{American Naturalist, 166} , 339-353.

Gunderson, A. R., \& Stillman, J. H. (2015). Plasticity in thermal
tolerance has limited potential to buffer ectotherms from global
warming. \emph{Proceedings of the Royal Society B: Biological Sciences,
282} (1808), 20150401. doi:10.1098/rspb.2015.0401

Harrison PM, Gutowsky LF, Martins EG, Patterson DA, Cooke SJ, Power M.
(2016). Temporal plasticity in thermal-habitat selection of burbot Lota
lota a diel-migrating winter-specialist. J Fish Biol. 2016
Jun;88(6):2111-29. doi: 10.1111/jfb.12990. Epub 2016 Apr 29. PMID:
27125426.

Kelly, M. W., Sanford, E., \& Grosberg, R. K. (2012). Limited potential
for adaptation to climate change in a broadly distributed marine
crustacean. \emph{Proceedings of the Royal Society B: Biological
Sciences, 279} (1727), 349-356. doi:doi:10.1098/rspb.2011.0542

Kurylyk, B. L.,~MacQuarrie, K. T. B.,~Linnansaari, T.,~Cunjak, R. A.,
and~Curry, R. A.~(2015)~Preserving, augmenting, and creating cold-water
thermal refugia in rivers: concepts derived from research on the
Miramichi River, New Brunswick (Canada).~\emph{Ecohydrol.}
,~8:~1095--~1108. doi:~10.1002/eco.1566.

Lande, R. (2009). Adaptation to an extraordinary environment by
evolution of phenotypic plasticity and genetic
assimilation.\emph{Journal of Evolutionary Biology, 22} (7), 1435-1446.

Lande, R. (2014). Evolution of phenotypic plasticity and environmental
tolerance of a labile quantitative character in a fluctuating
environment. \emph{Journal of Evolutionary Biology, 27} (5), 866-875.

MacLean, H. J., Sørensen, J. G., Kristensen, T. N., Loeschcke, V.,
Beedholm, K., Kellermann, V., \& Overgaard, J. (2019). Evolution and
plasticity of thermal performance: an analysis of variation in thermal
tolerance and fitness in 22 \emph{Drosophila}
species.\emph{Philosophical Transactions of the Royal Society B:
Biological Sciences, 374} (1778), 20180548.
doi:doi:10.1098/rstb.2018.0548

Padfield, D. \& Matheson, G. (2020). nls.multstart: robust non-linear
regression using AIC scores. R package version 1.2.0.
https://CRAN.R-project.org/package=nls.multstart

Padilla, D. K., \& Adolph, S. C. (1996). Plastic inducible morphologies
are not always adaptive: The importance of time delays in a stochastic
environment. \emph{Evolutionary Ecology, 10} (1), 105-117.
doi:10.1007/BF01239351

Pereira, R. J., Sasaki, M. C., \& Burton, R. S. (2017). Adaptation to a
latitudinal thermal gradient within a widespread copepod species: the
contributions of genetic divergence and phenotypic
plasticity.\emph{Proceedings of the Royal Society B: Biological
Sciences, 284} (1853), 20170236. doi:doi:10.1098/rspb.2017.0236

Phillips, B. L., Muñoz, M. M., Hatcher, A., Macdonald, S. L., Llewelyn,
J., Lucy, V., \& Moritz, C. (2016). Heat hardening in a tropical lizard:
geographic variation explained by the predictability and variance in
environmental temperatures. \emph{Functional Ecology, 30} (7),
1161-1168. doi:10.1111/1365-2435.12609

Pinheiro J, Bates D, DebRoy S, Sarkar D, R Core Team (2022).~\emph{nlme:
Linear and Nonlinear Mixed Effects Models} . R package version
3.1-155,~https://CRAN.R-project.org/package=nlme.

R Core Team (2021). R: A language and environment for statistical
computing. R Foundation for Statistical Computing, Vienna, Austria. URL
https://www.R-project.org/.

Rohr, J.R., Civitello, D.J., Cohen, J.M., Roznik, E.A., Sinervo, B. and
Dell, A.I. (2018), The complex drivers of thermal acclimation and
breadth in ectotherms. Ecol Lett, 21:
1425-1439.~https://doi.org/10.1111/ele.13107

Sgro, C. M., Overgaard, J., Kristensen, T. N., Mitchelll, K. A.,
Cockerell, F. E., \& Hoffmann, A. A. (2010). A comprehensive assessment
of geographic variation in heat tolerance and hardening capacity in
populations of \emph{Drosophila melanogaster} from eastern
Australia.\emph{Journal of Evolutionary Biology, 23} (11), 2484-2493.
doi:10.1111/j.1420-9101.2010.02110.x

Siljestam, M., \& Östman, Ö. (2017). The combined effects of temporal
autocorrelation and the costs of plasticity on the evolution of
plasticity. \emph{Journal of Evolutionary Biology, 30} (7), 1361-1371.
doi:10.1111/jeb.13114

Sultan, S. E., \& Spencer, H. G. (2002). Metapopulation structure favors
plasticity over local adaptation.\selectlanguage{english} . \emph{The American Naturalist, 160}
(2), 271-283. doi:10.1086/341015

Van Buskirk, J., \& Steiner, U. K. (2009). The fitness costs of
developmental canalization and plasticity. \emph{Journal of Evolutionary
Biology, 22} (4), 852-860.
doi:https://doi.org/10.1111/j.1420-9101.2009.01685.x

van Heerwaarden, B., Lee, R. F. H., Overgaard, J., \& Sgro, C. M.
(2014). No patterns in thermal plasticity along a latitudinal gradient
in \emph{Drosophila simulans} from eastern Australia. \emph{Journal of
Evolutionary Biology, 27} (11), 2541-2553. doi:10.1111/jeb.1251022

van Heerwaarden, B., Kellermann, V., \& Sgro, C. M. (2016). Limited
scope for plasticity to increase upper thermal limits. \emph{Functional
Ecology, 30} (12), 1947-1956. doi:10.1111/1365-2435.12687\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image1/image1}

\end{center}
\end{figure}

\textbf{Fig. 1.} Observed rates of plasticity (\selectlanguage{greek}\emph{λ\textsubscript{E}}\selectlanguage{english}
) in temperature tolerance for different classes of ectotherms animals.
Estimates obtained from a linear mixed effect model that controls for
the random effect of species as well as acclimation temperature and the
slope of the estimated exponential decay function at the final
measurement are given in Table 2.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image2/image2}

\end{center}
\end{figure}

\textbf{Fig. 2.} Observed rates of plasticity (\selectlanguage{greek}\emph{λ\textsubscript{E}}\selectlanguage{english}
) in temperature tolerance as a function of acclimation temperature for
different classes of ectothermic animals. Regression lines with a common
slope are based on the best fitting model (Table 2), and are given for
the median slope of the estimated exponential decay function at the
final acclimation time point, together with 95\% CI.

\textbf{Table 1.} 1 AICc comparisons of the three top candidate models
explaining variation in rates of plasticity (\selectlanguage{greek}\emph{λ\textsubscript{E}}\selectlanguage{english} )
in temperature tolerance among different classes of ectothermic animals.
'Slope' is the slope of the estimated exponential decay function at the
final measurement. Species identity is included as random intercepts in
all models.\selectlanguage{english}
\begin{longtable}[]{@{}lllll@{}}
\toprule
& K & AIC\textsubscript{C} & \selectlanguage{greek}Δ\selectlanguage{english}AIC\textsubscript{C} &
w\textsubscript{i}\tabularnewline
\midrule
\endhead
Acclimation temperature + Taxonomic class + Slope & 9 & -1031.8 & 0.00 &
0.73\tabularnewline
Acclimation temperature + Taxonomic class + Slope + Body mass & 10 &
-1029.8 & 2.00 & 0.27\tabularnewline
Class + Slope & 8 & -1019.1 & 12.70 & 0.00\tabularnewline
\bottomrule
\end{longtable}

\textbf{Table 2.} Summary of the best fitting linear mixed-effect model
(Table 1, re-fitted using REML) describing variation in rates of
plasticity (\selectlanguage{greek}\emph{λ\textsubscript{E}}\selectlanguage{english} ) in thermal tolerance among
different classes of ectothermic animals.\selectlanguage{english}
\begin{longtable}[]{@{}lll@{}}
\toprule
& Estimate & SE\tabularnewline
\midrule
\endhead
Fixed effects & &\tabularnewline
Amphibians & 0.053 & 0.005\tabularnewline
Reptiles & 0.032 & 0.008\tabularnewline
Insects & 0.011 & 0.007\tabularnewline
Fishes & 0.005 & 0.006\tabularnewline
Crustaceans & 0.004 & 0.009\tabularnewline
Acclimation temperature & 0.0006 & 0.0002\tabularnewline
Slope at end & 7.51 & 1.38\tabularnewline
Random effects (SD) & &\tabularnewline
Species & 0.017 &\tabularnewline
\bottomrule
\end{longtable}

\selectlanguage{english}
\FloatBarrier
\end{document}

