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 :