\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{Adaptive and non-adaptive plasticity in changing environments:
implications for sexual species with different life history strategies}



\author[1]{Daniel Romero-Mujalli}%
\author[1]{Markus Rochow}%
\author[1]{Sandra Kahl}%
\author[1]{Sofia Paraskevopoulou}%
\author[2]{Remco Folkertsma}%
\author[1]{Florian Jeltsch}%
\author[1]{Ralph Tiedemann}%
\affil[1]{University of Potsdam}%
\affil[2]{University of Potsdam Faculty of Mathematics and Natural Sciences}%


\vspace{-1em}



  
  \date{June 14, 2024}


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





\selectlanguage{english}
\begin{abstract}
Populations adapt to novel environmental conditions by genetic changes
or phenotypic plasticity. Plastic responses are generally faster and can
buffer fitness losses under variable conditions. Plasticity is typically
modelled as random noise and linear reaction norms that assume simple
one-to-one genotype-phenotype maps and no limits to the phenotypic
response. Most studies on plasticity have focused on its effect on
population viability. However, it is not clear, whether the advantage of
plasticity depends solely on environmental fluctuations or also on the
genetic and demographic properties (life histories) of populations. Here
we present an individual-based model and study the relative importance
of adaptive and non-adaptive plasticity for populations of sexual
species with different life histories experiencing directional
stochastic climate change. Environmental fluctuations were simulated
using differentially autocorrelated climatic stochasticity or noise
color, and scenarios of directional climate change. Non-adaptive
plasticity was simulated as a random environmental effect on trait
development, while adaptive plasticity as a linear, logistic, or
sinusoidal reaction norm. The last two imposed limits to the plastic
response and emphasized flexible interactions of the genotype with the
environment. Interestingly, this assumption led to (i) smaller
phenotypic than genotypic variance in the population and the coexistence
of polymorphisms, (ii) many-to-one genotype-phenotype map, and (iii) the
maintenance of higher genetic variation -- compared to linear reaction
norms and genetic determinism -- even when the population was exposed to
a constant environment for several generations. Limits to plasticity led
to genetic accommodation, when costs were negligible, and to the
appearance of cryptic variation when limits were exceeded. We found that
adaptive plasticity promoted population persistence under red noise
stochasticity and was particularly important for life histories with low
fecundity. Populations producing more offspring could cope with
environmental fluctuations solely by genetic changes or random
plasticity, unless environmental change was too fast.%
\end{abstract}%



\sloppy


Abstract

Populations adapt to novel environmental conditions by genetic changes
or phenotypic plasticity. Plastic responses are generally faster and can
buffer fitness losses under variable conditions. Plasticity is typically
modelled as random noise and linear reaction norms that assume simple
one-to-one genotype-phenotype maps and no limits to the phenotypic
response. Most studies on plasticity have focused on its effect on
population viability. However, it is not clear, whether the advantage of
plasticity depends solely on environmental fluctuations or also on the
genetic and demographic properties (life histories) of populations. Here
we present an individual-based model and study the relative importance
of adaptive and non-adaptive plasticity for populations of sexual
species with different life histories experiencing directional
stochastic climate change. Environmental fluctuations were simulated
using differentially autocorrelated climatic stochasticity or noise
color, and scenarios of directional climate change. Non-adaptive
plasticity was simulated as a random environmental effect on trait
development, while adaptive plasticity as a linear, logistic, or
sinusoidal reaction norm. The last two imposed limits to the plastic
response and emphasized flexible interactions of the genotype with the
environment. Interestingly, this assumption led to (\emph{i} ) smaller
phenotypic than genotypic variance in the population and the coexistence
of polymorphisms, (\emph{ii} ) many-to-one genotype-phenotype map, and
(\emph{iii} ) the maintenance of higher genetic variation -- compared to
linear reaction norms and genetic determinism -- even when the
population was exposed to a constant environment for several
generations. Limits to plasticity led to genetic accommodation, when
costs were negligible, and to the appearance of cryptic variation when
limits were exceeded. We found that adaptive plasticity promoted
population persistence under red noise stochasticity and was
particularly important for life histories with low fecundity.
Populations producing more offspring could cope with environmental
fluctuations solely by genetic changes or random plasticity, unless
environmental change was too fast.

Keywords

Individual-based models; phenotypic plasticity; reaction norms; limits;
many-to-one genotype-phenotype map; developmental canalization;
evolvability; noise color; environmental change; stochastic fluctuations

Introduction

A prevailing challenge in ecology and evolutionary biology is to
understand and predict species' responses to environmental change, such
as climate change (Chevin et al. 2010; Gonzalez et al. 2013).
Populations of different species are expected to respond to these
changes by adaptation to novel local environmental conditions, shifts in
their distributional range while tracking their preferred niche, or
local extinction (Franks et al. 2014; Wiens 2016). Particularly, when
movement opportunities are constrained, populations are expected to
either go extinct or to cope with novel conditions through adaptation by
genetic changes, or phenotypic plasticity. Phenotypic plasticity is
defined as the property of organisms sharing the same genotype to
produce different phenotypes, often in response to the local environment
(Pigliucci 2005; Reusch 2014; Sommer 2020). It constitutes the
non-heritable part of phenotypic variation and includes acclimation,
developmental plasticity, behavioral flexibility, learning, maternal
effects, epigenetics, and random noise (West-Eberhard 2003).

Among other factors, evolutionary rescue depends considerably on
demographic properties and the generation time of a species (Chevin et
al. 2010; Bell 2013). Thus, organisms that do not reproduce often or
produce relatively few progeny are expected to be vulnerable to
extinction under environmental change due to their lower supply of
beneficial mutations (Botero et al. 2015). Life history strategies where
a genetic response is limited, are expected to buffer fitness loss in a
novel environment through plastic responses, which are generally faster
(Lande 2009; Botero et al. 2015).Using individual-based models,
different life history strategies have been compared for their ability
to adapt to the local environment through genetic changes (e.g.,
Bj\selectlanguage{ngerman}örklund et al. 2009). In contrast, the role of plasticity --
particularly, the role of different types of phenotypic plasticity --
for persistence under environmental change, has not yet been thoroughly
investigated relative to different life history strategies. For example,
in species with similar generation time, those with relatively high
fecundity (``large clutch-size'') may rely less on plasticity as
compared to species with clutch-size limited to few offspring. Though it
was not the focus of their work, Björklund et al. (2009) observed in
their model that, all else equal, r-like life history strategists
persisted environmental change the longest (as compared to other life
history strategies) under scenarios of low heritability in which most
variability of the phenotypic trait was developed randomly (random
plasticity).

Phenotypic plasticity has long been considered important for organisms
experiencing fluctuating environmental conditions (Scheiner 1993; Via et
al. 1995). Yet, the modeling of phenotypic plasticity has not been a
straightforward task, since some features found in empirical research
remain elusive to current approaches (e.g., limits to plasticity, Murren
et al. 2015; many-to-one genotype-phenotype maps, Wagner 2005). Models
on the evolution of phenotypic plasticity suggest that plasticity
evolves in variable environmental conditions, when cues are reliable and
when costs are relatively low (Ghalambor et al. 2007; Chevin et al.
2010; Reed et al. 2010; Lande 2014; Ashander et al. 2016; Ergon and
Ergon 2016; Hendry 2016). However, year-to-year fluctuations of climatic
variables such as temperature can differ with regard to their serial
autocorrelation between consecutive time units (typically years; Laakso
et al. 2001, 2004; Schwager et al. 2006; Björklund et al. 2009). How
such environmental stochasticity may promote the degree and mode of
plastic responses has received less attention. Ecological models have
already shown that environmental stochasticity reduces long-term
population growth and that the type (i.e., the color) of the stochastic
noise differently affects population extinction risk. For instance,
Mustin et al. (2013) found that extinction risk is expected to be high
for populations experiencing directional climate change and inhabiting
climates with reddish (i.e., positively autocorrelated) stochasticity.
However, they did not consider scenarios of negatively autocorrelated
stochastic noise (blue noise), nor the effect of plasticity on
population persistence.

Here, we present an individual-based eco-evolutionary model to study the
relative importance of adaptive and non-adaptive plasticity for
populations of sexual species with different life histories,
experiencing scenarios of various rates of directional climate change
and different types of environmental stochasticity (noise color). This
paper deals with organisms experiencing rapid environmental change and
hence does not consider the evolution of plasticity itself. Instead, it
focuses on so far understudied aspects of plasticity by (\emph{i} )
specifically testing for the influence of limits to plasticity,
(\emph{ii} ) exploring developmental flexibility ; and (\emph{iii} )
evaluating the performance of different types of plasticity and life
histories under environmental change (e.g., we ask: what type of
plasticity is preferred under a particular life history?). Assuming
limits to the plastic response may add realism to the modeling of
plasticity. For example, when costs of plasticity are negligible, as it
seems to be the case for the majority of organisms (Murren et al. 2015),
a model assuming perfect sensing and no limits may arrive at perfectly
adapted phenotypes under all circumstances, while in nature, the plastic
response is limited to a definite range of environmental conditions,
beyond which perfect plasticity is physiologically impossible
(Wiesenthal et al. 2018). This is not the case, however, for linear
reaction norms, where the perfect plasticity is prevented only when
assuming costs to plasticity or biased perception of environmental cues.
Most theoretical work in the literature has focused on understanding the
adaptive nature of phenotypic plasticity (Via and Lande 1985; Nussey et
al. 2007; Chevin et al. 2010; Lande 2014), mainly using linear reaction
norms and simple one-to-one genotype-phenotype mapping (Chevin et al.
2010; Reed et al. 2010; Lande 2014). However, plasticity can result --
at least for quantitative traits -- from a complex relationship between
genotype and phenotype, with the developmental system responding
flexibly to internal (genotype) and external inputs (environment)
(Laland 2015). To this end, we compare linear reaction norms with
alternative plasticity types, including a flexible developmental system.
As a consequence, multiple genotypes can have the same phenotype and are
mutationally interconnected (many-to-one genotype-phenotype map, Wagner
2008; Ahnert 2017; Aguilar-Rodríguez et al. 2018). This assumption leads
to interesting evolutionary and ecological outcomes that, to our
knowledge, are not captured by previous models, and that may be more in
line with what empirical evidence suggests.

\emph{}

2 Methods

To study the effect of adaptive and non-adaptive phenotypic plasticity
on population persistence under scenarios of environmental change, we
developed an eco-evolutionary individual-based model (IBM) of a
geographically isolated panmictic population of a sexual species with
non-overlapping generations experiencing stochastic directional climate
change. The focus was on studying the ability of a population to adapt
to its local environment (no migration was possible). This modeling
setup could resemble a fish population inhabiting a lake, or a plant or
animal population inhabiting a highly fragmented environment where
movement opportunities are constrained. Populations could differ in
fecundity and intrinsic population dynamics (different life history
strategies). The model also allows for different forms of variation in
environmental stochasticity or noise color: uncorrelated white noise
typical for terrestrial locations; positively autocorrelated red noise,
which had been found in coastal and marine habitats, Vasseur and Yodzis
2004; and negatively autocorrelated blue noise. Blue noise is less
common, but recent evaluations of climate spectral exponents suggest
that temperature has turned bluer (i.e., tends towards more negatively
autocorrelated stochasticity) in most continents in the last century
(García-Carreras and Reuman 2011).

The model was created using the freely available software platform
Netlogo 6.0.2 (Wilensky 1999) and is available for download from
(place-link-here). A full description of the model that follows the ODD
(Overview, Design, concepts, and Details) protocol (Grimm et al. 2006,
2010) can be found in Appendix A. Below, only model features that were
used in this study are explained. The sequence of model operations was:
set initial environment and population (assumed to be locally adapted);
update phenotypic response, check for degree of adaptation (as fitness
proxy), computation of fecundity, reproduction of adults, inheritance,
die-off of adults, check for extinction and update of the environmental
state before repeating the loop (Fig. 1).

\textbf{Fig. 1} Flowchart diagram of the sequence of model operations.

2.1 Landscape

The environment imposed a phenotypic optimum
\selectlanguage{greek}\emph{θ\textsubscript{t}}\selectlanguage{english}(hereafter, environmental optimum) which could
change at constant speed every generation depending on the simulated
scenario of environmental change. Thus, \selectlanguage{greek}\emph{θ\textsubscript{t} =
θ\textsubscript{0} + η t}\selectlanguage{english}determined the directional trend of the optimum
\selectlanguage{greek}\emph{θ\textsubscript{t}}\selectlanguage{english} in a deterministic environment (no
stochasticity). The parameter\selectlanguage{greek}\emph{θ\textsubscript{0}}\selectlanguage{english} was the initial
environmental optimum (when\emph{t = 0} ) and \selectlanguage{greek}\emph{η} \selectlanguage{english}was the rate of
environmental change. By varying the parameter \selectlanguage{greek}\emph{η} \selectlanguage{english}we simulated
different scenarios of directional climate change (e.g., no change,
slow, medium, rapid climate change). Stochastic colored noise around
\selectlanguage{greek}\emph{θ\textsubscript{t}}\selectlanguage{english} was implemented to simulate different
scenarios of environmental stochasticity (Fig. 2). This method has been
recommended for the simulation of directional climate change scenarios
(Kopp and Matuszewski 2014; Vincenzi 2014).

\textbf{Fig. 2} Example of scenarios of directional climate change and
environmental stochasticity as simulated in the model. These scenarios
resemble fluctuations of climatic variables, such as temperature, as
measured for different kinds of habitats. Different forms of stochastic
fluctuations (noise color) of the environmental optimum
(\selectlanguage{greek}\emph{θ\textsubscript{t}}\selectlanguage{english} ) were simulated. They differ in their
temporal autocorrelation, i.e., no autocorrelation (white), negative
autocorrelation (blue), or positive autocorrelation (red). The dashed
lines illustrate different rates of directional climate change (no
change \selectlanguage{greek}\emph{η = 0}\selectlanguage{english} , slow to rapid change, scenarios of \selectlanguage{greek}\emph{η
\textgreater{} 0}\selectlanguage{english} ). Different scenarios of \selectlanguage{greek}\emph{η} \selectlanguage{english}were explored per
scenario of stochastic colored noise (environmental autocorrelation).
The color bar shows the range of values explored for the level of
autocorrelation (\selectlanguage{greek}\emph{α} , \selectlanguage{english}see methods).

Stochasticity according to colored noise was implemented such that the
environmental optimum was redefined as \selectlanguage{greek}\emph{θ\textsubscript{t} =
θ*\textsubscript{t} + ϕ\textsubscript{t}}\selectlanguage{english}
where\selectlanguage{greek}\emph{θ*\textsubscript{t}}\selectlanguage{english} gave the directional trend of the mean
environmental optimum as specified above and \selectlanguage{greek}\emph{ϕ\textsubscript{t}=
αϕ\textsubscript{t-1} + βξ\textsubscript{t}}\selectlanguage{english} the stochastic noise. The
parameter \selectlanguage{greek}\emph{α} \selectlanguage{english}governed the level of environmental autocorrelation
and therefore allowed for different forms of stochasticity or noise
color as in Bj\selectlanguage{ngerman}örklund et al (2009): \selectlanguage{greek}\emph{-1 \textless{} α \textless{}
0}\selectlanguage{english} , blue noise; \selectlanguage{greek}\emph{α = 0}\selectlanguage{english} , white noise, and\selectlanguage{greek}\emph{0 \textless{} α
\textless{} 1}\selectlanguage{english} , red noise (Fig. 2). Several scenarios of noise color
(values of \selectlanguage{greek}\emph{α} \selectlanguage{english}) were explored, ranging from negatively
autocorrelated environmental conditions or blue noise over uncorrelated
(white noise) to positively autocorrelated environmental conditions or
red noise (see Table I). The parameter\selectlanguage{greek}\emph{β = σ} \selectlanguage{english}\(\sqrt{1-\alpha^{2}}\)
was the adjusted environmental variance for all degrees of
autocorrelation, as in (Schwager et al. 2006), and
\selectlanguage{greek}\emph{σ\textsuperscript{2}}\selectlanguage{english} \emph{= 1} was the environmental variance.
The parameter \selectlanguage{greek}\emph{ξ\textsubscript{t}}\selectlanguage{english} was a random value, normally
distributed with zero mean and unity as variance.

2.2 The population

Individuals in the population were characterized by sex, stage (whether
adult or juvenile), degree of adaptation (fitness proxy, depends on the
ecological phenotype and the environment), fecundity (influenced by the
degree of adaptation and the population density), and an ecological
phenotype (evolving trait). The ecological phenotype \emph{z} had its
genetic component \emph{a} defined by \emph{L} diploid loci with
additive effects, and its environmental component \emph{e} determining
the contribution of phenotypic plasticity in the development of the
trait. Thus, \emph{z\textsubscript{i} = a\textsubscript{i} + e} , was
the ecological phenotype of individual \emph{i} , where
\emph{a\textsubscript{i}}was the sum of allelic values at the adaptation
locus. We assumed one diploid locus affecting the phenotypic trait. At
the beginning of the simulation, the population was composed of \emph{N}
individuals at carrying capacity (\emph{N = K = 1000} ) and locally
adapted. This means that for each individual organism, alleles coding
for its ecological phenotype were drawn from a normal distribution with
mean equal to the environmental optimum\selectlanguage{greek}\emph{θ\textsubscript{0}}\selectlanguage{english} and
variance \emph{V = VG / 2} , where \emph{VG = 0.2} was the initial
genetic variance present in the population.

2.3 Phenotypic plasticity

Because we wanted to compare traditional approaches with alternative
ones imposing limits to plasticity (objective 1, see Introduction), four
different types or scenarios of phenotypic plasticity were implemented
(Fig. 3): random noise, linear reaction norm, sinusoidal reaction norm,
and logistic reaction norm. In the model, phenotypic plasticity affected
the environmental component \emph{e} of the ecological phenotype. Random
noise has been the most common method in eco-evolutionary IBMs for the
representation of an environmental effect on the development of the
phenotypic trait (Romero-Mujalli et al. 2019b). In this model, we
consider random noise as part of non-adaptive phenotypic plasticity.
Thus, we could compare the effect of non-adaptive and adaptive
phenotypic plasticity on population persistence. In the model, random
plasticity assumed \emph{e} to be random and normally distributed with
zero mean and variance \selectlanguage{greek}\emph{VE = σ\textsuperscript{2}}\selectlanguage{english} , the
environmental variance.

\textbf{Fig. 3} Forms of non-adaptive (random) and adaptive (linear,
sinusoidal, logistic) phenotypic plasticity implemented in the model.
The black line indicates the moving phenotypic
optimum\selectlanguage{greek}\emph{θ\textsubscript{t}}\selectlanguage{english} as given by the environment. Empty
circles show the phenotypic response of the organism.

For the modeling of adaptive phenotypic plasticity, we considered three
approaches: the linear reaction norm without a limit, which is the most
common approach in the literature (Chevin et al. 2010; Lande 2014;
Romero-Mujalli et al. 2019b); and two other approaches that account for
limits to phenotypic plasticity (Fig. 3A). In the literature, most
research has focused on costs of plasticity, and very little on its
limits (Murren et al. 2015).

In our model, the linear reaction norm defined the environmental
component as \selectlanguage{greek}\emph{e = bθ\textsubscript{t}}\selectlanguage{english} ; where \emph{b} is the
slope or degree of plasticity, and \selectlanguage{greek}\emph{θ\textsubscript{t}}\selectlanguage{english} the
environmental optimum at time \emph{t} (in generations). Furthermore,
the linear reaction norm was assumed to be shallower (\emph{b = 0.5} )
and in the direction of the phenotypic optimum \selectlanguage{greek}\emph{θ\textsubscript{t}}\selectlanguage{english}
. Hence, the linear plastic response was adaptive, but that there was a
misfit arising from an imperfect development of the trait (Lande 2009;
Chevin et al. 2010; Ashander et al. 2016). It is important to mention
that perfect sensing of the environment was assumed for all scenarios of
adaptive phenotypic plasticity. Reliability of sensing, i.e., the
correlation between sensed environmental cues and environmental
variables affecting fitness, has been investigated by Reed et al.
(2010), Ashander et al. (2016) and Ergon and Ergon (2016), but was
beyond the scope of our work. This was decided in favor of highlighting
fundamental differences between traditional (linear reaction norms) and
alternative (logistic and sinusoidal) approaches.

The two other methods, logistic and sinusoidal, were designed based on
observations from stress tolerance responses for some physiological
traits (Jordan and Deaton 1999; Ara\selectlanguage{ngerman}újo et al. 2014, 2016; Solan and
Whiteley 2016; Wiesenthal et al. 2018), and on behavioral traits. Their
plastic response was assumed to be developmentally flexible (Laland et
al. 2015), relying on feedback with the environment, and in the
direction of the environmental optimum (adaptive). This means that they
allowed for stable functioning of the phenotype (close to the optimum)
despite of the variation at the genetic level. Furthermore, in contrast
to linear reaction norms, they enable multiple genotypes to have the
same phenotype, which is a widespread feature of genotype-phenotype maps
(Ahnert 2017; Aguilar-Rodríguez et al. 2018; Payne and Wagner 2019).
These methods differed from each other only in the condition that
determined the phenotype produced when the limit is exceeded, and were
implemented as follows:

\selectlanguage{greek}\emph{e = ΩΔE}\selectlanguage{english} ; where \selectlanguage{greek}\emph{Ω} \selectlanguage{english}was always positive and shaped the
plastic response. It was given by \selectlanguage{greek}\emph{Ω = sin(}\selectlanguage{english} |\selectlanguage{greek}\emph{ΔE}\selectlanguage{english} |\emph{)} .
The term \selectlanguage{greek}\emph{ΔE}\selectlanguage{english} indicated the amount of change with respect to the
reference environment \selectlanguage{greek}\emph{θ\textsuperscript{R} = a}\selectlanguage{english} , such
that\selectlanguage{greek}\emph{ΔE = θ\textsubscript{t} -- θ\textsuperscript{R}}\selectlanguage{english} . The
parameter\emph{a} was the genetic component of the phenotype.

For sinusoidal phenotypic plasticity, if the argument of the sine
function was greater than \selectlanguage{greek}\emph{π} , \selectlanguage{english}the environmental component
\emph{e}was set to \emph{0} (the organism fails to develop a plastic
response). An example could be snails subject to salinity stress. If the
change is too large (compared to the reference environment where
plasticity is not needed), snails fail to develop enough physiological
response to counter and balance osmotic pressure (e.g., Jordan and
Deaton 1999; Wiesenthal et al. 2018).

For logistic phenotypic plasticity, if the argument of the sine function
was greater than \selectlanguage{greek}\emph{π / 2}\selectlanguage{english} ; the term \selectlanguage{greek}\emph{ΔE}\selectlanguage{english} was set to \selectlanguage{greek}\emph{ΔE =
π / 2}\selectlanguage{english} such that a maximum response was reached (saturation). This could
resemble plant species exposed to different light conditions. After some
point of increasing light intensity, a maximum thickness will be
reached, and the plant's leaves would not grow any thicker (Wilson and
Cooper 1969).

According to the sequence of model operations (Figure 1), the plastic
response always followed the environmental change. This implementation
implies phenotypic traits being flexible enough such that adult
individuals could still respond to environmental change. However, some
characters cannot be further modified after a sensitive period early in
life (constant characters, Lande 2014; Romero-Mujalli et al. 2019b),
such that adults experiencing a new environment will no longer be able
to adjust their matured phenotype. In order to evaluate the impact of
the ontogenetic phase of plasticity (flexible or labile plasticity vs.
fixed during ontogeny or constant characters after an early sensitive
period) on population persistence, we also study the effect of
phenotypic plasticity when the plastic response precedes the
environmental change.

2.4 Degree of adaptation

After developing the phenotype, adult individuals in the population were
subject to selection according to the following Gaussian function
(Burger and Lynch 1995):

\begin{equation}
w_{i}=e^{\frac{{-(z_{i}-\ \theta_{t})}^{2}}{2\gamma^{2}}}\nonumber \\
 \end{equation}

where \emph{w\textsubscript{i}} was the degree of adaptation (fitness
proxy) of individual \emph{i} , and \selectlanguage{greek}\emph{γ} , \selectlanguage{english}the width of the function
(strength of selection). The closer the ecological
phenotype\emph{z\textsubscript{i}} was to the optimum
\selectlanguage{greek}\emph{θ\textsubscript{t}}\selectlanguage{english} , the better the individual coped with the
environmental conditions.

2.5 Fecundity and life history strategies

The fecundity of individuals in the population was scaled according to
their degree of adaptation after considering density dependence effects,
and follows a Ricker-type model as in Bj\selectlanguage{ngerman}örklund et al (2009):

\begin{equation}
{w^{\prime}}_{i}={w_{i}e}^{\psi(1-\frac{N}{K})}\nonumber \\
 \end{equation}

where \emph{w'\textsubscript{i}} was the fecundity of individual\emph{i}
, \emph{N} was the population size and \emph{K} , the carrying capacity.
The parameter \selectlanguage{greek}\emph{ψ} \selectlanguage{english}described the strength of the density dependence
effect (also referred here as density compensation). The higher \selectlanguage{greek}\emph{ψ}
, \selectlanguage{english}the stronger was the effect of increased population density \emph{N} .
Here we implemented three levels of \selectlanguage{greek}\emph{ψ} \selectlanguage{english}as in Bj\selectlanguage{ngerman}örklund et al.
(2009): \emph{0.5} , \emph{1.8} and \emph{2.5} . These three values
produced fundamentally different population dynamics and therefore cover
a wide range of different outcomes due to \selectlanguage{greek}\emph{ψ} \selectlanguage{english}(Fig. 4A and B).

\textbf{Fig. 4} Life history strategies as implemented in the model. (A)
expected fecundity (\selectlanguage{greek}\emph{λ} \selectlanguage{english}) per reproductive couple for different
values of population size (carrying capacity \emph{K = 1000} ). When the
population size is low, resources are plenty, and well adapted couples
can contribute their best in number of offspring (fecundity) for the
next generation. (B) Population dynamics (carrying capacity \emph{K =
1000} ) per life history strategy (values of \selectlanguage{greek}\emph{ψ} \selectlanguage{english}) in a static
environment (no climate change, no stochasticity).

2.6 Reproduction

Adult individuals mated randomly with others of opposite sex, with
replacement for males only (i.e., lottery polygyny, males could
participate in more than one reproductive event). The fecundity of the
reproductive couple \selectlanguage{greek}\emph{λ} \selectlanguage{english}was equal to the sum of the scaled fitness
of the partners \emph{i} ,\emph{j} :

\selectlanguage{greek}\emph{λ = w'\textsubscript{i} + w'\textsubscript{j}}\selectlanguage{english}

Each couple produced a number of offspring randomly drawn from a Poisson
distribution with expectancy \selectlanguage{greek}\emph{λ} .

\selectlanguage{english}2.7 Inheritance

After reproduction, each offspring inherited one strain copy or
haplotype from each parent. Only the genetic component (or breeding
value) \emph{a} was inherited. Allele values of the inherited haplotype
were picked randomly from the corresponding parental locus. Each
haplotype mutated with probability \selectlanguage{greek}\emph{μ} \selectlanguage{english}of mutation per locus. In
case that a mutation occurred, its effect was randomly drawn from a
normal distribution with zero mean and variance equal to the effect size
of mutations \emph{MV} , which was an input parameter (Vincenzi 2014).
The assumption of a Gaussian distribution is consistent with analysis of
mutation effects (Lynch and Walsh 1998; Martin et al. 2006).

All adults died after reproduction (non-overlapping generations), and
the environmental state was updated before repeating the loop.

2.8 Simulation experiments

With the model we compared different implementations of plasticity
(objective 1) and studied the effect of these types of plasticity
(non-adaptive and adaptive phenotypic plasticity) on the persistence of
populations with different life histories under scenarios of directional
stochastic environmental change (objective 2). Probability of
persistence was computed as the proportion of times a persistence event
(or success) occurs from the total number of observations (the 100
replicates). Error bars were estimated from the binomial distribution. A
total of 100 replicates lasting 250 generations each were performed for
different combinations of rate of directional climate change, degree of
temporal autocorrelation (noise color), type of phenotypic plasticity
(including genetic determinism), and density dependence effect (life
history strategy) (see table I). To illustrate selected parameter
values: under no stochasticity, \selectlanguage{greek}\emph{η = 0.02}\selectlanguage{english} imposed a directional
trend on the optimum phenotype of 0.02 phenotypic units per generation
(or per year, assuming one generation a year). This means that after 250
generations of simulation time, the optimum
\selectlanguage{greek}\emph{θ\textsubscript{t}}\selectlanguage{english}would have changed in five phenotypic units,
and a population with no density compensation (i.e., \selectlanguage{greek}\emph{ψ} \selectlanguage{english}\emph{= 1}
) could avoid extinction via evolutionary rescue or some sort of
phenotypic plasticity. To our experimental setup, the contribution of
mutations was important. Otherwise, populations failed to cope with the
changing environmental conditions (Fig. S1, Appendix B: Supplementary
material). It is important to highlight that, in our model, limits of
plasticity were reached well before the simulation ends. Following the
same example, an initially well adapted individual
(\selectlanguage{greek}\emph{z\textsubscript{i} = θ\textsubscript{0}}\selectlanguage{english} ) reached its limits of
plasticity (i.e., saturation, by logistic plasticity; decay in ability
of plastic response, by sinusoidal plasticity) at
\selectlanguage{greek}\emph{θ\textsubscript{t}}\selectlanguage{english} \emph{= 1.58} . This corresponds to generation
\emph{t = 80} .

To complement the findings on life history strategies of intermediate
and strong density dependence effects, additional scenarios of rate of
environmental change and of stochastic fluctuations in carrying capacity
were implemented to unravel whether, in the model, life history
strategies differently tolerate fluctuations in K. Specifically, each
time unit of model iteration (or generation), the value of carrying
capacity was drawn from a normal distribution centered in \emph{K} and
variance \emph{sd} \textsuperscript{2}, where \emph{sd} was a percentage
of \emph{K} . For instance, if \emph{K = 1000} , the scenario of 10\%
of\emph{K} means \emph{sd = 100} . This is, around 70\% of the cases,
the value of carrying capacity is expected to occur within the interval
(\emph{K -- sd} , \emph{K + sd} ).

Simulation experiments were also evaluated under scenarios of lower
mutation rate, and mutational effects according to the model of slightly
deleterious mutations (Ohta 1973; Eyre-Walker et al. 2002;
Romero-Mujalli et al. 2019a), which imposed higher genetic constraints.
The analysis of data and plotting was performed in r v3.5.2.

\textbf{Table I} Parameter values and description\selectlanguage{english}
\begin{longtable}[]{@{}lll@{}}
\toprule
\textbf{Parameter} & \textbf{Value} &
\textbf{Description}\tabularnewline
\midrule
\endhead
\emph{K} & 1000 & Carrying capacity\tabularnewline
\selectlanguage{greek}\emph{γ} & 2.2 & Strength of selection. Moderate selection (Bj\selectlanguage{ngerman}örklund et
al. 2009)\tabularnewline
\selectlanguage{greek}\emph{σ\textsuperscript{2}}\selectlanguage{english} & 1 & Variance of the stochastic
environment\tabularnewline
\selectlanguage{greek}\emph{θ\textsubscript{0}}\selectlanguage{english} & 0 & Initial environmental
optimum\tabularnewline
\selectlanguage{greek}\emph{η} & 0, 0.01, 0.02, 0.025, 0.03, 0.04 & Rate of environmental
change\tabularnewline
\selectlanguage{greek}\emph{α} & -0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8 & Level of
environmental autocorrelation (scenarios of
stochasticity)\tabularnewline
\selectlanguage{greek}\emph{ψ} & 0.5, 1.0, 1.8, 2.5 & Density dependence effect (life history
strategy) (as in Bj\selectlanguage{ngerman}örklund et al. 2009)\tabularnewline
\emph{z} & evolving trait & Ecological phenotype\tabularnewline
\emph{b} & 0.5 & Slope of the linear reaction norm\tabularnewline
\emph{VG} & 0.2 & Initial genetic variance present in the population
(Vincenzi 2014)\tabularnewline
\emph{L} & 1 & Number of loci per evolving trait\tabularnewline
\selectlanguage{greek}\emph{μ} & 10\textsuperscript{-3}, 10\textsuperscript{-4} & Mutation
rate*\tabularnewline
\emph{MV} & 0.2 & Variance of the distribution of mutations fitness
effect size (Vincenzi 2014)\tabularnewline
\emph{t} & 100, 250 & Time limit per simulation, in
generations\tabularnewline
\emph{sd} & 10 and 20\% the value of \emph{K} & Standard deviation of
the stochastic carrying capacity\tabularnewline
\bottomrule
\end{longtable}

* The value of the mutation rate was picked according to results of a
simulation model on the evolution of the mutation rate after 300
generations of a population experiencing stochastic environmental
conditions (no climate change, and no plasticity, Romero-Mujalli et al.
2019a). In addition this value is within the range of mutation rates
used in other simulation models (reviewed in Romero-Mujalli et al.
2019b).

3. Results

3.1 Implications of assuming limits and flexible development for the
modeling of phenotypic plasticity

Considering the degree of adaptation under directional deterministic
environmental change (no stochastic noise, typical of laboratory
experiments), any form of adaptive phenotypic plasticity was of
advantage over non-adaptive plasticity and absence of plasticity (i.e.,
genetic determinism) (Fig. 5A). When the slope equals one, linear
reaction norms yield perfect plasticity, while the ``perfect match'' of
the phenotype with the environmental optimum depends on whether the
limits are exceeded or not for those forms of plasticity that account
for limits (Fig. 5A).

\textbf{Fig. 5} (A) Mean degree of adaptation (fitness proxy) in
populations with and without adaptive phenotypic plasticity experiencing
directional deterministic environmental change (rate of change \selectlanguage{greek}\emph{η =
0.03}\selectlanguage{english} , density dependence effect \selectlanguage{greek}\emph{ψ = 1.0}\selectlanguage{english} , mutation rate \selectlanguage{greek}\emph{µ
= 0.001}\selectlanguage{english} ). A scenario of no adaptation (\selectlanguage{greek}\emph{µ = 0}\selectlanguage{english} ) was also
simulated, for comparison. The simulation was run for 100 generations.
(B) Genetic variance (average \selectlanguage{ngerman}± standard error from 100 replicates)
present in a population experiencing constant environmental conditions
(no climate change, no stochasticity; for example, laboratory
conditions) per scenario of phenotypic plasticity. (C) Time series of
the genotypic and phenotypic variances present in a population
experiencing deterministic directional environmental change (no
stochasticity), under the scenarios of adaptive sinusoidal and logistic
phenotypic plasticity (upper panel). When the population is locally
adapted (initial part), there is less phenotypic than genetic variation
(\emph{developmental canalization} , at the population level, Posadas
and Carthew 2014). Under the assumption of adaptive sinusoidal
plasticity, as the environment changes and the population is pushed
towards its limits of plasticity, cryptic phenotypic variation arises,
and the mean fitness (here, degree of adaptation) of the population
reduces (lower panel). In linear reaction norms, the genotypic and
phenotypic variances are directly proportional.

Interestingly, adaptive sinusoidal and logistic plasticity led to stable
functioning phenotypes (close to the optimum) despite of variation at
the genetic level (\emph{developmental canalization} at the population
level, Posadas and Carthew 2014) (Figure 5B, C). A consequence of this
property is that they maintained higher genetic variance than the
alternative methods even when the population was exposed to a constant
environment (e.g., laboratory conditions) for several generations (Fig.
5B). Additionally, in the model assuming a sinusoidal reaction norm,
cryptic variation appeared in the population when it was pushed towards
the limits (Fig. 5C). Cryptic variation refers to genetic variation that
normally has little or no effect on phenotypic variation but that under
atypical conditions generates phenotypic variation (Paaby and Rockman
2014).

3.2 Phenotypic plasticity: life histories under environmental change

Under scenarios of weak density dependence effect (\selectlanguage{greek}\emph{ψ = 0.5}\selectlanguage{english} ), in
which breeding females in the population could produce relatively few
offspring, adaptive phenotypic plasticity played a major role promoting
persistence as compared to organisms with higher \selectlanguage{greek}\emph{ψ} \selectlanguage{english}(Fig. 6, 7).
Particularly for this life history strategy under slow climate change,
adaptive phenotypic plasticity became of high importance promoting
adaptation under positively autocorrelated environmental stochasticity
(red noise). When the rate of change \selectlanguage{greek}\emph{η} \selectlanguage{english}was too rapid, adaptive
phenotypic plasticity (particularly, linear reaction norm and logistic
plasticity) became of advantage for all simulated environmental
conditions of noise color (Fig. 6).

\textbf{Fig. 6} The effect of non-adaptive (random) and adaptive
(linear, sinusoidal, and logistic) phenotypic plasticity on probability
of persistence (100 replicates, 250 generations each) of a population in
which breeding females are limited to produce relatively few offspring
(\selectlanguage{greek}\emph{ψ = 0.5}\selectlanguage{english} , weak density dependence effect). The linear reaction
norm had a slope \emph{b = 0.5} . A scenario of genetic determinism
(narrow sense heritability \emph{h\textsuperscript{2} = 1} ) was also
simulated, for comparison. The color bar illustrates the color of the
stochastic noise.

\textbf{Fig. 7} Effect of non-adaptive (random) and adaptive
(sinusoidal, logistic, and linear) phenotypic plasticity on the
probability of persistence (100 replicates, 250 generations each) of a
population of intermediate and strong density dependence effect (\selectlanguage{greek}\emph{ψ
= 1.8 and 2.5}\selectlanguage{english} , respectively). The linear reaction norm had a
slope\emph{b = 0.5} . A scenario of genetic determinism (narrow sense
heritability \emph{h\textsuperscript{2} = 1} ) was also simulated. The
color bar illustrates the color of the stochastic noise. Scenarios
of\selectlanguage{greek}\emph{η \textless{} 0.03}\selectlanguage{english} were not shown, because all treatments led
to maximum probability of persistence.

In contrast, results of strong and very strong density compensation
(\selectlanguage{greek}\emph{ψ = 1.8 and 2.5}\selectlanguage{english} , respectively) always persisted under
conditions of only environmental stochasticity (no directional climate
change) and relatively slow to medium rate of environmental change
(\selectlanguage{greek}\emph{η = 0, 0.01, 0.02}\selectlanguage{english} ), regardless of the type of environmental
stochasticity (noise color). For these life history strategies, genetic
determinism and all forms of plasticity performed equally well (Fig. 7).

Under scenarios of relatively rapid directional climate change, adaptive
linear and logistic phenotypic plasticity performed the best for
organisms of intermediate and high density dependence effects (\selectlanguage{greek}\emph{ψ =
1.8}\selectlanguage{english} and \selectlanguage{greek}\emph{ψ = 2.5}\selectlanguage{english} ) (Fig. 8, \selectlanguage{greek}\emph{η = 0.04}\selectlanguage{english} ). The life history
strategy with weak density compensation was not included in this
analysis, since its populations always went extinct, except for the
scenarios of linear and logistic plasticity. On the other hand, adaptive
sinusoidal, non-adaptive random phenotypic plasticity, and genetic
determinism showed similar performance across all scenarios of rapid
rate of directional climate change (Fig. 7).

The performance of random plasticity was notably prominent for organisms
of intermediate and high density dependence effects under scenarios of
lower mutation rates (\selectlanguage{greek}\emph{μ} \selectlanguage{english}= 10\textsuperscript{-4}, Fig. 8).
Similar results were observed when considering \selectlanguage{greek}\emph{μ} \selectlanguage{english}=
10\textsuperscript{-3} and mutational effects according to the model of
slightly deleterious mutations (Ohta 1973; Eyre-Walker et al. 2002), as
in Romero-Mujalli et al. (2019b) (Fig. S2, Appendix B: Supplementary
material). When genetic mechanisms producing novel variation are somehow
constrained, these organisms can cope with changing environmental
conditions relaying on non-adaptive phenotypic plasticity only (Fig. 8).
This ability was not observed for organisms of weak density
compensation, typical of large mammals and some bird species. However,
populations of organisms with strong density compensation (\selectlanguage{greek}\emph{ψ =
2.5}\selectlanguage{english} ) were considerably more vulnerable to extinction due to stochastic
fluctuations in the carrying capacity (Fig. 9).

\textbf{Fig. 8} Relative importance of forms of non-adaptive (random)
and adaptive (sinusoidal, logistic and linear) phenotypic plasticity
affecting persistence of populations differing in levels of density
compensation (\selectlanguage{greek}\emph{ψ} \selectlanguage{english}) and experiencing moderate rate of directional
stochastic environmental change (\selectlanguage{greek}\emph{η = 0.01}\selectlanguage{english} , \selectlanguage{greek}\emph{ψ = 0.5}\selectlanguage{english} ;
\selectlanguage{greek}\emph{η = 0.025}\selectlanguage{english} , \selectlanguage{greek}\emph{ψ = 1.8}\selectlanguage{english} and \selectlanguage{greek}\emph{ψ = 2.5;}\selectlanguage{english} 100 replicates,
250 generations each) under scenarios of low mutation rate (\selectlanguage{greek}\emph{μ} \selectlanguage{english}=
10\textsuperscript{-4}). A higher \selectlanguage{greek}\emph{η} , \selectlanguage{english}results for \selectlanguage{greek}\emph{ψ = 0.5}\selectlanguage{english}
were not shown, because scenarios of random and sinusoidal plasticity,
and genetic determinism (narrow sense heritability
\emph{h\textsuperscript{2} = 1} ) went usually extinct. The color bar
illustrates the color of the stochastic noise.

\textbf{Fig. 9} Probability of persistence (100 replicates, 250
generations each) per life history strategy (levels of \selectlanguage{greek}\emph{ψ} \selectlanguage{english}) under
stochastic fluctuations in both, the carrying capacity (\emph{sd =
10}and \emph{20\%} of \emph{K} ) and the phenotypic optimum (white
noise, no directional climate change).

When the plastic response precedes the environmental change (i.e.,
change after a sensitive period of plastic trait development), the
dynamic can considerably change for those life history strategies of
weak density compensation (\selectlanguage{greek}\emph{ψ = 0.5}\selectlanguage{english} ) under blue stochastic noise
(Fig. S3, Appendix B: Supplementary material). Whether plastic responses
precede or follow the environmental change seems to have little impact
under life history strategies with stronger density compensation (Fig.
S3, Appendix B: Supplementary material).

4. Discussion

The objectives of this study were: (\emph{i} ) to evaluate the effect of
assuming plasticity as a developmentally flexible phenotypic response
with limits, and (\emph{ii} ) to assess the relative importance of
adaptive and non-adaptive plasticity for populations of sexual species
with different life history strategy experiencing scenarios of
directional climate change and environmental stochasticity (noise
color). The simulated environmental conditions, though simplified as in
every model, mimic realistic expected scenarios of environmental climate
change (Bj\selectlanguage{ngerman}örklund et al. 2009; Kopp and Matuszewski 2014; Vincenzi
2014). Our results show that the relative importance of phenotypic
plasticity varies among life history strategies. Furthermore, they show
that the advantage of non-adaptive and adaptive forms of phenotypic
plasticity on population persistence depends on the type of
environmental stochasticity (noise color) and the rate of directional
climate change. In addition, the advantage of adaptive and non-adaptive
forms of plasticity depends on whether plasticity precedes or follows
the environmental change. Finally, assuming limits to plasticity lead to
genetic accommodation (in the absence of costs to plasticity) and to the
appearance of cryptic genetic variation when limits are exceeded.
Moreover, assuming plasticity as the result of genotypes responding
flexibly to feedback from their environment leads to developmental
canalization at the population level (Posadas and Carthew 2014),
many-to-one genotype-phenotype map (Wagner 2008; Aguilar-Rodríguez et
al. 2018) promoting the coexistence of polymorphism, and higher
maintained genetic variation, potentially increasing evolvability, as
compared to traditional approaches (e.g., linear reaction norms).

4.1 On the limits and the model implementations of adaptive phenotypic
plasticity

Our work explicitly accounts for limits to plasticity. Though a thorough
analysis of limits to plasticity was beyond the scope of this study,
there are some aspects worth consideration. Linear reaction norms have
no limits and their evolution can theoretically result in organisms
displaying perfectly adapted phenotypes under all environmental
contexts, if no cost is imposed, which is unrealistic. Plastic responses
certainly have energetic demands (i.e., costs), but also limits to
plasticity -- in the absence of costs -- can result from the properties
and characteristics of the underlying machinery producing the phenotypic
response. Therefore, linear adaptive phenotypic plasticity can
overestimate the probability of persistence. In our simulations this was
prevented by having a slope of 0.5 instead of 1. Even so, linear
adaptive phenotypic plasticity had the best performance (in terms of
increasing population persistence) in comparison to the other methods of
adaptive plasticity. In the model, it took approximately 80, 60 and 40
generations or years to cross the limits of plasticity for a rate of
directional change of 0.02, 0.03 and 0.04, respectively. Thus, for the
population to persist the 250 generations under such scenarios of
environmental change, genetic changes need to necessarily follow to
enable survival beyond plasticity limits. Only linear adaptive
phenotypic plasticity, as simulated in the model, could -- for some life
history strategies, \selectlanguage{greek}\emph{ψ = 1.8}\selectlanguage{english} and \emph{2.5} -- lead to persistence
of the population with minimal, if any, genetic changes. Thus, the
existence of limits to plasticity (without the consideration of costs;
i.e., when costs are negligible, see Murren et al. 2015) can lead to a
transition from an ``environmentally induced'' to a ``genetically
encoded'' state of a trait (genetic accommodation, Sommer 2020). The
process of genetic accommodation occurs when phenotypic variants that
are environmentally induced, become genetically determined by natural
selection (West-Eberhard 2003; Schlichting and Wund 2014), and has been
observed in nature (Schlichting and Wund 2014; Kulkarni et al. 2017).
Moreover, our model shows that when the environmental change pushes
towards the limits of plasticity, cryptic phenotypic variation may
arise. This means higher phenotypic than genotypic variance, a ratio
that could be measured in nature. This feature comes along with a
reduction in mean fitness of the population (in the model, degree of
adaptation) and could potentially be used as an early warning signal for
the inability of a population to sustain environmental change (Boettiger
et al. 2013).

In our model, forms of phenotypic plasticity with limits are assumed to
enable diverse functional responses through flexible interactions of the
genotype with the environment (Laland et al. 2015; van Gestel and
Weissing 2016). This assumption leads to, first, stable functioning of
the phenotype (close to the optimum) despite the variation at the
genetic level (\emph{developmental canalization} at the population
level, Posadas and Carthew 2014). Second, multiple genotypes can have
the same phenotype (many-to-one genotype-phenotype map, Ahnert 2017;
Aguilar-Rodr\selectlanguage{ngerman}íguez et al. 2018), which does not occur under linear
reaction-norms, thus enabling coexistence of polymorphism when different
genotypes have the same fitness. Third, genotypes with the same
phenotype are mutationally interconnected, such that small mutations can
transform these genotypes into one another without altering the
phenotype (robustness), or lead to new phenotypes (evolvability), at the
same time (Wagner 2008; Ahnert 2017; Aguilar-Rodríguez et al. 2018;
Payne and Wagner 2019). Moreover, it can act as a mechanism that
maintains genetic variation, potentially increasing evolvability (van
Gestel and Weissing 2016). In our implementation of adaptive plasticity
this is evidenced by the sinusoidal and logistic plasticity scenarios
resulting in the highest maintained genetic variation -- higher than
linear reaction norms and genetic determinism -- even when the
population is exposed to a constant environment (e.g., laboratory
conditions) for several generations. This agrees with observations on
the maintenance of adaptive capacity in populations (often clonal
populations) kept in enclosed environments (e.g., chemostats) for
several generations (Maharjan et al. 2006; Fussmann et al. 2007).
Modelling approaches have been limited by the tendency of the trait
variance to unrealistically decline to zero over time (Merico et al.
2014; Romero-Mujalli et al. 2019b). Our study shows that a simulation of
trait dynamics more in line with developmental systems responding
flexibly to internal (genotype) and external inputs (environment)
enabling highly diverse functional responses, as done for logistic and
sinusoidal plasticity, could help to overcome this limitation. So far,
our implementation offers a phenomenological alternative only. Future
theoretical work should focus on unraveling the mechanisms that make the
above-discussed phenotypic response (trait dynamics) possible.
Furthermore, and in contrast to other forms of plasticity, flexible
development of the phenotype can provide an adaptive response in a novel
environment that need not to have been pre-screened by earlier
selection, which is particularly true for learning (Laland et al. 2015;
Watson and Szathmáry 2016; Romero-Mujalli et al. 2017).

4.2 The relative importance of adaptive and non-adaptive phenotypic
plasticity differs between life history strategies

According to our model, any form of adaptive phenotypic plasticity is of
advantage over non-adaptive plasticity under controlled systematic
changes of environmental conditions, as it is the case of laboratory
experiments. However, under stochastic environmental fluctuations
typical of natural habitats, adaptive phenotypic plasticity is
particularly of advantage when facing positively autocorrelated
environmental fluctuations (i.e., red noise) and slow directional
environmental change; and for a broader spectrum of climatic
fluctuations when there is a faster trend in the mean environmental
optimum. Particularly, all forms of adaptive phenotypic plasticity were
superior to lack of plasticity for life history strategies with
relatively slow growth rate (weak density compensation). This life
history strategy resembles characteristics of many vertebrate species
(e.g., birds and mammals). In contrast, for life history strategies with
strong density dependence effects and therefore, fast growth rate,
genetic determinism and all forms of plasticity performed equally well.
For these life history strategies, only under relatively high rates of
directional climate change, some forms of adaptive phenotypic plasticity
(i.e., linear and logistic phenotypic plasticity) and, interestingly,
random non-adaptive plasticity, were of advantage. The difference in the
relative importance of adaptive and non-adaptive plasticity among life
history strategies was even more evident under scenarios of genetic
constraints (in our model, scenarios of relatively low mutation rate,
and of slightly deleterious mutations fitness effects). Under such
scenarios, life history strategies with strong density compensation
could cope with the changing environmental conditions relying on random
non-adaptive plasticity only. This was not the case for life history
strategies with weak density compensation, which strongly depended on
adaptive plasticity for their adaptation to the local environment.
Therefore, if organisms with different life history strategies, as
simulated in this model, would experience equivalent environmental
fluctuations and rates of directional climate change, those where
breeding females are limited to few offspring are expected to experience
stronger selection for the development of mechanisms of adaptive
phenotypic plasticity, and stronger selection to expand the limits of
their plasticity mechanisms. Populations of organisms with faster growth
rate can rely on non-adaptive plasticity (translating into bet-hedging,
Donaldson-Matasci et al. 2013) for a broader range of environmental
fluctuations. As linear and logistic phenotypic plasticity does not
generally outperform non-adaptive plasticity, these adaptive plastic
responses are not expected to evolve for life history strategies
producing relatively large numbers of offspring, unless they experience
relatively rapid rates of directional environmental change. A further
life history parameter promoting the evolution of adaptive phenotypic
plasticity is longevity/generation time, leading to a limited genetic
response (Forsman 2014, not tested in our study).

The importance of adaptive phenotypic plasticity for organisms
experiencing directional change of the mean environmental optimum, as
inferred in our simulations, may equally apply to dispersing and sessile
organisms. Dispersing organisms may experience gradual changes in the
mean environmental optimum and will benefit from developing adaptive
forms of phenotypic plasticity as they expand their range, especially if
density dependence is weak, as it occurs in mammal and bird species.
Similarly, sessile organisms exposed to seasonal changes in the mean
environmental optimum will also benefit from adaptive forms of
phenotypic plasticity, particularly if their longevity is high (Borges
2008). Examples are plant species inhabiting temperate regions
(Chmielewski and Rötzer 2001), as well as plants experiencing regular
yearly cycles of rain, drought and fire at the equator (Fajardo et al.
2005).

Despite of their good performance across scenarios of directional
stochastic environmental change (as tested in this study), populations
with strong density dependence are susceptible to extinction under
stochastic fluctuations in environmental quality (i.e., stochastic
carrying capacity). Fluctuations in the carrying capacity can result
from fire, resource contamination, or human land use, among other
factors (Anderson et al. 2015). In contrast, populations with weak
density compensation are less impacted by fluctuations in the carrying
capacity of their environment. This would, for example, suggest that
insect populations are more impacted by fluctuations in habitat quality
(e.g., because of land use) than mammals. It remains however to be
investigated, in how far this observation depends on the specific model
implementation.

4.3 Phenotypic plasticity and environmental noise: when adaptive
plasticity hinders evolution

To date, most studies of phenotypic plasticity (and its evolution) under
stochastic environmental fluctuations have been focused on the level of
correlation among the environmental optimum and the environmental cues
sensed by organisms (Reed et al. 2010; Ashander et al. 2016; Ergon and
Ergon 2016). They show that adaptive phenotypic plasticity can only
evolve under positive correlation of cues with the environmental optimum
(environmental predictability). It is important to note that this type
of predictability is not the same as the predictable (red noise)
year-to-year pattern of our simulations. Specifically, this study adds
that adaptive forms of phenotypic plasticity can decrease extinction
risk under positively autocorrelated environmental stochasticity (red
noise), while extinction risk is expected to be high in the absence of
plasticity (Mustin et al. 2013). However, if the environment changes
after the sensitive period when the trait is susceptible to be modified
by plasticity, those traits initially developed by adaptive plasticity
can lead to a disadvantage under negatively autocorrelated environmental
stochasticity (blue noise), which agrees with results from evolution
experiments (Hallsson and Björklund 2012). This applies for constant
characters, i.e., phenotypic traits plastic early in ontogeny, resulting
in mature phenotypes which cannot be further modified by the environment
(Lande 2014; Romero-Mujalli et al. 2019b). In our model, such
disadvantage of adaptive plasticity under blue noise occurred only for
life history strategies with weak density compensation. Whether the
plastic response precedes or follows the environmental change seems to
have little impact under life history strategies with stronger density
compensation. Thus, for life history strategies with weak density
compensation under blue noise, adaptive phenotypic plasticity amplifies
phenotype-environment mismatches due to the negative autocorrelation in
the environmental stochasticity and, therefore, hinders evolution in the
long run, whilst non-adaptive phenotypic plasticity and genetic
determinism increase population persistence under such scenarios.
Therefore, life history strategies with weak density compensation
relying on any form of adaptive phenotypic plasticity under blue
stochastic environmental conditions are expected to be more vulnerable
-- than those with stronger density compensation -- unless they are able
to adjust their plasticity strategy. However, adaptive plasticity
(particularly, linear and logistic phenotypic plasticity) becomes of
advantage for all conditions of stochastic noise -- even when the
environmental change occurs after the sensitive period -- when the rate
of environmental change is relatively rapid. In our model, rapid rates
of directional change of 0.03 and 0.04 change the environment beyond the
range captured by the stochasticity (considering a two-sigma-effect,
95\%) of the original reference environment after 100 and 75 generations
(or years), respectively. A ``fast'' rate of environmental change is not
necessarily an absolute measurement. It should be related to
characteristics -- life history strategy -- of the population, as shown
in this study, and to the strength of selection (Burger and Lynch 1995).

5 Conclusions

Adaptive phenotypic plasticity promotes population persistence under
positively autocorrelated (red noise) environmental stochasticity and
slow climate change, and for a broader range of fluctuations when the
rate of directional change is faster. This form of plasticity is
particularly important for life history strategies in which breeding
females have a limited number of offspring (low fecundity, and hence
slow population growth rate, typical of many vertebrate species: birds
and mammals). Organisms producing more offspring per female may cope
with environmental fluctuations relying only on genetic changes or
random plasticity, unless the rate of environmental change is relatively
high. Whenever plasticity precedes the environmental change, typical of
constant characters, those life history strategies with weak density
compensation will experience high risk of extinction in bluish habitats,
unless they can adjust their plasticity into a bet-hedging strategy.

Models employing linear reaction norms may overestimate persistence, if
they do not consider limits of plasticity. Furthermore, limits to
plasticity lead to genetic accommodation when costs are negligible and
to the exposure of cryptic genetic variation when the plastic response
is pushed towards the limits of plasticity. In addition, the modelling
of plasticity as developmental systems relying on genotypes interacting
flexibly with their environment enables coexistence of polymorphisms and
highly diverse functional responses, leading to high maintained genetic
variation, many-to-one genotype-phenotype map, and to developmental
canalization at the population level.

In this work, the mechanisms that shape the limits of adaptive
plasticity were not explicitly modeled. Empirical work is needed to
unravel molecular mechanisms that may dictate the limits of plastic
responses.

References

Aguilar-Rodríguez, J., L. Peel, M. Stella, A. Wagner, and J. L. Payne.
2018. The architecture of an empirical genotype-phenotype map. Evolution
72:1242--1260. John Wiley \& Sons, Ltd.

Ahnert, S. E. 2017. Structural properties of genotype-phenotype maps. J.
R. Soc. Interface 14:20170275. The Royal Society.

Anderson, C., Z. Jovanoski, I. Towers, and H. Sidhu. 2015. A simple
population model with a stochastic carrying capacity.

Araújo, C. V. M., E. N. V. Rodríguez, D. Salvatierra, L. A.
Cedeño-Macias, V. C. Vera-Vera, M. Moreira-Santos, and R. Ribeiro. 2016.
Attractiveness of food and avoidance from contamination as conflicting
stimuli to habitat selection by fish. Chemosphere 163:177--183.

Araújo, C. V. M., C. Shinn, M. Moreira-Santos, I. Lopes, E. L. G.
Espíndola, and R. Ribeiro. 2014. Copper-driven avoidance and mortality
in temperate and tropical tadpoles. Aquat. Toxicol. 146:70--75.

Ashander, J., L.-M. Chevin, and M. L. Baskett. 2016. Predicting
evolutionary rescue via evolving plasticity in stochastic environments.
Proc. R. Soc. B Biol. Sci. 283:20161690.

Bell, G. 2013. Evolutionary rescue and the limits of adaptation. Philos.
Trans. R. Soc. B Biol. Sci. 368.

Björklund, M., E. Ranta, V. Kaitala, L. A. Bach, P. Lundberg, and N.
Chr. Stenseth. 2009. Quantitative Trait Evolution and Environmental
Change. PLoS ONE 4:e4521.

Boettiger, C., N. Ross, and A. Hastings. 2013. Early warning signals:
the charted and uncharted territories. Theor. Ecol. 6:255--264.

Borges, R. M. 2008. Plasticity comparisons between plants and animals:
Concepts and mechanisms. Plant Signal. Behav. 3:367--375.

Botero, C. A., F. J. Weissing, J. Wright, and D. R. Rubenstein. 2015.
Evolutionary tipping points in the capacity to adapt to environmental
change. Proc. Natl. Acad. Sci. 112:184.

Burger, R., and M. Lynch. 1995. Evolution and Extinction in a Changing
Environment: A Quantitative-Genetic Analysis. Evolution 49:151--163.

Chevin, L.-M., R. Lande, and G. M. Mace. 2010. Adaptation, Plasticity,
and Extinction in a Changing Environment: Towards a Predictive Theory.
PLoS Biol 8:e1000357.

Chmielewski, F.-M., and T. Rötzer. 2001. Response of tree phenology to
climate change across Europe. Agric. For. Meteorol. 108:101--112.

Donaldson-Matasci, M. C., C. T. Bergstrom, and M. Lachmann. 2013. When
Unreliable Cues Are Good Enough. Am. Nat. 182:313--327.

Ergon, T., and R. Ergon. 2016. When three traits make a line: evolution
of phenotypic plasticity and genetic assimilation through linear
reaction norms in stochastic environments. J. Evol. Biol. 30:486--500.

Eyre-Walker, A., P. D. Keightley, N. G. C. Smith, and D. Gaffney. 2002.
Quantifying the Slightly Deleterious Mutation Model of Molecular
Evolution. Mol. Biol. Evol. 19:2142--2149.

Fajardo, L., V. González, J. M. Nassar, P. Lacabana, C. A. Portillo Q.,
F. Carrasquel, and J. P. Rodríguez. 2005. Tropical Dry Forests of
Venezuela: Characterization and Current Conservation Status. Biotropica
37:531--546.

Forsman, A. 2014. Rethinking phenotypic plasticity and its consequences
for individuals, populations and species. Heredity 115:276.

Franks, S. J., J. J. Weber, and S. N. Aitken. 2014. Evolutionary and
plastic responses to climate change in terrestrial plant populations.
Evol. Appl. 7:123--139.

Fussmann, G. F., M. Loreau, and P. A. Abrams. 2007. Eco-evolutionary
dynamics of communities and ecosystems. Funct. Ecol. 21:465--477.

García-Carreras, B., and D. C. Reuman. 2011. An empirical link between
the spectral colour of climate and the spectral colour of field
populations in the context of climate change. J. Anim. Ecol.
80:1042--1048.

Ghalambor, C. K., J. K. Mckay, S. P. Carroll, and D. N. Reznick. 2007.
Adaptive versus non-adaptive phenotypic plasticity and the potential for
contemporary adaptation in new environments. Funct. Ecol. 21:394--407.

Gonzalez, A., O. Ronce, R. Ferriere, and M. E. Hochberg. 2013.
Evolutionary rescue: an emerging focus at the intersection between
ecology and evolution. Philos. Trans. R. Soc. B Biol. Sci. 368:20120404.

Grimm, V., U. Berger, F. Bastiansen, S. Eliassen, V. Ginot, J. Giske, J.
Goss-Custard, T. Grand, S. K. Heinz, G. Huse, A. Huth, J. U. Jepsen, C.
Jørgensen, W. M. Mooij, B. Müller, G. Pe'er, C. Piou, S. F. Railsback,
A. M. Robbins, M. M. Robbins, E. Rossmanith, N. Rüger, E. Strand, S.
Souissi, R. A. Stillman, R. Vabø, U. Visser, and D. L. DeAngelis. 2006.
A standard protocol for describing individual-based and agent-based
models. Ecol. Model. 198:115--126.

Grimm, V., U. Berger, D. L. DeAngelis, J. G. Polhill, J. Giske, and S.
F. Railsback. 2010. The ODD protocol: A review and first update. Ecol.
Model. 221:2760--2768.

Hallsson, L. R., and M. Björklund. 2012. Selection in a fluctuating
environment leads to decreased genetic variation and facilitates the
evolution of phenotypic plasticity. J. Evol. Biol. 25:1275--1290. John
Wiley \& Sons, Ltd.

Hendry, A. P. 2016. Key Questions on the Role of Phenotypic Plasticity
in Eco-Evolutionary Dynamics. J. Hered. 107:25--41.

Jordan, P. J., and L. E. Deaton. 1999. Osmotic regulation and salinity
tolerance in the freshwater snail Pomacea bridgesi and the freshwater
clam Lampsilis teres. Comp. Biochem. Physiol. A. Mol. Integr. Physiol.
122:199--205.

Kopp, M., and S. Matuszewski. 2014. Rapid evolution of quantitative
traits: theoretical perspectives. Evol. Appl. 7:169--191.

Kulkarni, S. S., R. J. Denver, I. Gomez-Mestre, and D. R. Buchholz.
2017. Genetic accommodation via modified endocrine signalling explains
phenotypic divergence among spadefoot toad species. Nat. Commun. 8:993.

Laakso, J., V. Kaitala, and E. Ranta. 2001. How does environmental
variation translate into biological processes? Oikos 92:119--122. John
Wiley \& Sons, Ltd.

Laakso, J., V. Kaitala, and E. Ranta. 2004. Non-linear biological
responses to environmental noise affect population extinction risk.
Oikos 104:142--148. John Wiley \& Sons, Ltd.

Laland, K. N., T. Uller, M. W. Feldman, K. Sterelny, G. B. Müller, A.
Moczek, E. Jablonka, and J. Odling-Smee. 2015. The extended evolutionary
synthesis: its structure, assumptions and predictions. Proc. R. Soc. B
282:20151019.

Lande, R. 2009. Adaptation to an extraordinary environment by evolution
of phenotypic plasticity and genetic assimilation. J. Evol. Biol.
22:1435--1446.

Lande, R. 2014. Evolution of phenotypic plasticity and environmental
tolerance of a labile quantitative character in a fluctuating
environment. J. Evol. Biol. 27:866--875.

Lynch, M., and B. Walsh. 1998. Genetics and Analysis of Quantitative
Traits. Sinauer.

Maharjan, R., S. Seeto, L. Notley-McRobb, and T. Ferenci. 2006. Clonal
Adaptive Radiation in a Constant Environment. Science 313:514.

Martin, G., T. Lenormand, and C. Goodnight. 2006. The fitness effect of
mutations across environments: a survey in light of fitness landscape
models. Evolution 60:2413--2427.

Merico, A., G. Brandt, S. L. Smith, and M. Oliver. 2014. Sustaining
diversity in trait-based models of phytoplankton communities. Front Ecol
Evol 2:59.

Murren, C. J., J. R. Auld, H. Callahan, C. K. Ghalambor, C. A.
Handelsman, M. A. Heskel, J. G. Kingsolver, H. J. Maclean, J. Masel, H.
Maughan, D. W. Pfennig, R. A. Relyea, S. Seiter, E. Snell-Rood, U. K.
Steiner, and C. D. Schlichting. 2015. Constraints on the evolution of
phenotypic plasticity: limits and costs of phenotype and plasticity.
Heredity 115:293--301.

Mustin, K., C. Dytham, T. G. Benton, and J. M. J. Travis. 2013. Red
noise increases extinction risk during rapid climate change. Divers.
Distrib. 19:815--824.

Nussey, D. H., A. J. Wilson, and J. E. Brommer. 2007. The evolutionary
ecology of individual phenotypic plasticity in wild populations. J.
Evol. Biol. 20:831--844.

Ohta, T. 1973. Slightly Deleterious Mutant Substitutions in Evolution.
Nature 246:96.

Paaby, A. B., and M. V. Rockman. 2014. Cryptic genetic variation:
evolution's hidden substrate. Nat. Rev. Genet. 15:247.

Payne, J. L., and A. Wagner. 2019. The causes of evolvability and their
evolution. Nat. Rev. Genet. 20:24--38.

Pigliucci, M. 2005. Evolution of phenotypic plasticity: where are we
going now? Trends Ecol. Evol. 20:481--486.

Posadas, D. M., and R. W. Carthew. 2014. MicroRNAs and their roles in
developmental canalization. Curr. Opin. Genet. Dev. 27:1--6.

Reed, T. E., R. S. Waples, D. E. Schindler, J. J. Hard, and M. T.
Kinnison. 2010. Phenotypic plasticity and population viability: the
importance of environmental predictability. Proc. R. Soc. B Biol. Sci.,
doi: 10.1098/rspb.2010.0771.

Reusch, T. B. H. 2014. Climate change in the oceans: evolutionary versus
phenotypically plastic responses of marine animals and plants. Evol.
Appl. 7:104--122.

Romero-Mujalli, D., J. Cappelletto, E. A. Herrera, and Z. Tárano. 2017.
The effect of social learning in a small population facing environmental
change: an agent-based simulation. J. Ethol. 35:61--73.

Romero-Mujalli, D., F. Jeltsch, and R. Tiedemann. 2019a. Elevated
mutation rates are unlikely to evolve in sexual species, not even under
rapid environmental change. BMC Evol. Biol. 19:175.

Romero-Mujalli, D., F. Jeltsch, and R. Tiedemann. 2019b.
Individual-based modeling of eco-evolutionary dynamics: state of the art
and future directions. Reg. Environ. Change 1--12.

Scheiner, S. M. 1993. Genetics and Evolution of Phenotypic Plasticity.
Annu. Rev. Ecol. Syst. 24:35--68.

Schlichting, C. D., and M. A. Wund. 2014. Phenotypic plasticity and
epigenetic marking: an assessment of evidence for genetic accommodation.
Evolution 68:656--672.

Schwager, M., K. Johst, and F. Jeltsch. 2006. Does Red Noise Increase or
Decrease Extinction Risk? Single Extreme Events versus Series of
Unfavorable Conditions. Am. Nat. 167:879--888.

Solan, M., and N. Whiteley. 2016. Stressors in the Marine Environment:
Physiological and ecological responses; societal implications. Oxford
University Press, Oxford, UK.

Sommer, R. J. 2020. Phenotypic Plasticity: From Theory and Genetics to
Current and Future Challenges. Genetics 215:1.

van Buskirk, J., and U. K. Steiner. 2009. The fitness costs of
developmental canalization and plasticity. J. Evol. Biol. 22:852--860.

van Gestel, J., and F. J. Weissing. 2016. Regulatory mechanisms link
phenotypic plasticity to evolvability. Sci. Rep. 6:24524.

Vasseur, D. A., and P. Yodzis. 2004. The color of environmental noise.
Ecology 85:1146--1152.

Via, S., R. Gomulkiewicz, G. D. Jong, S. M. Scheiner, C. D. Schlichting,
and P. H. V. Tienderen. 1995. Adaptive phenotypic plasticity: consensus
and controversy. Trends Ecol. Evol. 10:212--217.

Via, S., and R. Lande. 1985. Genotype-Environment Interaction and the
Evolution of Phenotypic Plasticity. Evolution 39:505--522.

Vincenzi, S. 2014. Extinction risk and eco-evolutionary dynamics in a
variable environment with increasing frequency of extreme events. J. R.
Soc. Interface 11.

Wagner, A. 2008. Robustness and evolvability: a paradox resolved. Proc.
R. Soc. B Biol. Sci. 275:91--100. Royal Society.

Wagner, A. 2005. Robustness and Evolvability in Living Systems:
Princeton University Press.

Watson, R. A., and E. Szathmáry. 2016. How Can Evolution Learn? Trends
Ecol. Evol. 31:147--157.

West-Eberhard, M. J. 2003. Developmental Plasticity and Evolution.
Oxford University Press, New York.

Wiens, J. J. 2016. Climate-Related Local Extinctions Are Already
Widespread among Plant and Animal Species. PLOS Biol. 14:e2001104.

Wiesenthal, A. A., C. Müller, and J.-P. Hildebrandt. 2018. Potential
modes of range shifts in euryhaline snails from the Baltic Sea and fresh
water lakes in northern Germany. Hydrobiologia 811:339--350.

Wilensky, U. 1999. NETLOGO. Centre for Connected Learning and
Computer-Based Modelling, Northwestern University, Evanston, IL.

Wilson, D., and J. P. Cooper. 1969. Effect of Light Intensity During
Growth on Leaf Anatomy and Subsequent Light-Saturated Photosynthesis
Among Contrasting Lolium Genotypes. New Phytol. 68:1125--1135.

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

