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).
Fig. 1 Flowchart diagram of the sequence of model operations.
2.1 Landscape
The environment imposed a phenotypic optimum θt(hereafter, environmental optimum) which could change at constant speed
every generation depending on the simulated scenario of environmental
change. Thus, θt = θ0 + η tdetermined the directional trend of the optimum θt in a deterministic environment (no stochasticity). The parameterθ0 was the initial environmental optimum (whent = 0 ) and η was the rate of environmental change. By
varying the parameter η we simulated different scenarios of
directional climate change (e.g., no change, slow, medium, rapid climate
change). Stochastic colored noise around θt 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).
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
(θt ) 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 η = 0 , slow to rapid change, scenarios of η
> 0 ). Different scenarios of η were explored per
scenario of stochastic colored noise (environmental autocorrelation).
The color bar shows the range of values explored for the level of
autocorrelation (α , see methods).
Stochasticity according to colored noise was implemented such that the
environmental optimum was redefined as θt =
θ*t + ϕt whereθ*t gave the directional trend of the mean
environmental optimum as specified above and ϕt= αϕt-1 + βξt the stochastic
noise. The parameter α governed the level of environmental
autocorrelation and therefore allowed for different forms of
stochasticity or noise color as in Björklund et al (2009): -1
< α < 0 , blue noise; α = 0 , white noise, and0 < α < 1 , red noise (Fig. 2). Several
scenarios of noise color (values of α ) 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β = σ \(\sqrt{1-\alpha^{2}}\) was the adjusted environmental
variance for all degrees of autocorrelation, as in (Schwager et al.
2006), and σ2 = 1 was the environmental
variance. The parameter ξt 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 z had its
genetic component a defined by L diploid loci with
additive effects, and its environmental component e determining
the contribution of phenotypic plasticity in the development of the
trait. Thus, zi = ai + e , was the
ecological phenotype of individual i ,
where aiwas 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 N individuals at
carrying capacity (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θ0 and variance V = VG / 2 , where 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 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 e to be random and normally distributed with
zero mean and variance VE = σ2 , the
environmental variance.
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θt 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 e = bθt ; where b is the slope
or degree of plasticity, and θt the environmental
optimum at time t (in generations). Furthermore, the linear
reaction norm was assumed to be shallower (b = 0.5 ) and in the
direction of the phenotypic optimum θt . 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ú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:
e = ΩΔE ; where Ω was always positive and shaped the
plastic response. It was given by Ω = sin( │ΔE │) .
The term ΔE indicated the amount of change with respect to the
reference environment θR = a , such thatΔE = θt – θR . The parametera was the genetic component of the phenotype.
For sinusoidal phenotypic plasticity, if the argument of the sine
function was greater than π , the environmental component ewas set to 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 π / 2 ; the term ΔE was set to ΔE =
π / 2 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 wi was the degree of adaptation (fitness
proxy) of individual i , and γ , the width of the function
(strength of selection). The closer the ecological phenotypezi was to the optimum θt ,
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örklund et al (2009):
\begin{equation}
{w^{\prime}}_{i}={w_{i}e}^{\psi(1-\frac{N}{K})}\nonumber \\
\end{equation}where w’i was the fecundity of individuali , N was the population size and K , the carrying
capacity. The parameter ψ described the strength of the density
dependence effect (also referred here as density compensation). The
higher ψ , the stronger was the effect of increased population
density N . Here we implemented three levels of ψ as in
Björklund et al. (2009): 0.5 , 1.8 and 2.5 . These
three values produced fundamentally different population dynamics and
therefore cover a wide range of different outcomes due to ψ (Fig.
4A and B).
Fig. 4 Life history strategies as implemented in the model. (A)
expected fecundity (λ ) per reproductive couple for different
values of population size (carrying capacity 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 K =
1000 ) per life history strategy (values of ψ ) 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 λ was equal to the sum of the scaled fitness
of the partners i ,j :