Skip to main content
Open Access

Species’ Range Dynamics Affect the Evolution of Spatial Variation in Plasticity under Environmental Change

Abstract

While clines in environmental tolerance and phenotypic plasticity along a single species’ range have been reported repeatedly and are of special interest in the context of adaptation to environmental changes, we know little about their evolution. Recent empirical findings in ectotherms suggest that processes underlying dynamic species’ ranges can give rise to spatial differences in environmental tolerance and phenotypic plasticity within species. We used individual-based simulations to investigate how plasticity and tolerance evolve in the course of three scenarios of species’ range shifts and range expansions on environmental gradients. We found that regions of a species’ range that experienced a longer history or larger extent of environmental change generally exhibited increased plasticity or tolerance. Such regions may be at the trailing edge when a species is tracking its ecological niche in space (e.g., in a climate change scenario) or at the front edge when a species expands into a new habitat (e.g., in an expansion/invasion scenario). Elevated tolerance and plasticity in the distribution center was detected when asymmetric environmental change (e.g., polar amplification) led to a range expansion. However, tolerance and plasticity clines were transient and slowly flattened out after range dynamics because of genetic assimilation.

Online enhancements:   appendixes, supplemental table and figures. Dryad data: https://dx.doi.org/10.5061/dryad.2nc0bn1.

Introduction

Species exhibit remarkable abilities to survive in variable environments, which manifests as environmental tolerance—a phenomenon that caught the interest of biologists early on (e.g., Grinnell 1917; Elton 1927; Hutchinson 1957). A species’ environmental tolerance can be broadly defined as its ecological niche, an important conceptual tool to understand the geographical distribution of species (Guisan and Zimmermann 2000; Essl et al. 2009; Slatyer et al. 2013) and to predict their response to environmental changes (Hijmans and Graham 2006; Valladares et al. 2014). Environmental tolerance can be realized through phenotypic plasticity as the capacity to express multiple environment-dependent phenotypes (Chevin et al. 2010). Alternatively, tolerance can result from the ability of a phenotype to perform in multiple environments, such as when morphological or physiological characteristics provide access to a broader range of resources (e.g., Roughgarden 1972). The capacity of a genotype to produce adapted phenotypes is, however, rarely perfect, and organisms often have an environment in which they perform best (e.g., Eppley 1972; Huey and Kingsolver 1989). The relationship between environmental variation and organismal performance thus often results in a unimodal curve with performance decreasing away from an optimal environmental condition. Yet not only do species differ from each other in the breadth of their niche, they also differ in their environmental tolerance and plasticity among populations within their range (e.g., Macdonald and Chinnappa 1989; Woods et al. 2012; Bennett et al. 2015; Lancaster 2016; Toftegaard et al. 2016; Reger et al. 2018). For example, terrestrial ectotherms often exhibit clines in environmental tolerance with broader thermal tolerances at high latitudes compared with species or populations near the equator (Addo-Bediako et al. 2000; Gaston 2009; Sunday et al. 2012; Lancaster et al. 2015; Lancaster 2016). The meta-analysis of Sunday et al. (2011) revealed between-species differences in thermal tolerance for 239 terrestrial ectotherm species (including arthropods, amphibians, and reptiles), with high-latitude species having a broader thermal tolerance than low-latitude species. Lancaster (2016) additionally assembled within-species thermal tolerance data for 17 insect species demonstrating that populations located closer to the poles had larger thermal tolerances than populations at more equatorial positions. While inter- and intraspecies spatial variation in tolerance breadth and plasticity are well documented, a comprehensive understanding of the underlying biological processes has not yet been achieved.

Divergence in tolerance and niche breadth can result from evolutionary processes, given that environmental tolerance is a heritable trait with additive genetic variance for fitness (Via and Lande 1985; Lynch and Gabriel 1987; Chevin and Lande 2011). Selection pressures affecting tolerance evolution mainly stem from the environmental variability that genotypes experience, across space or over time. In temporally variable environments (Lynch and Gabriel 1987; Lande 2014) or in structured populations with gene flow between distinct habitats (Via and Lande 1985; Sultan and Spencer 2002; van Buskirk 2017), environmental tolerance is expected to increase when genotypes encounter multiple environments over time. In contrast, lower environmental variability lessens the fitness benefits of broad environmental tolerance and limits its evolution when a trade-off is present (Lynch and Gabriel 1987; Padilla and Adolph 1996; Reed et al. 2010; Ezard et al. 2014). Maintaining environmental tolerance via phenotypic plasticity may also entail physiological or metabolic costs, further opposing its evolution (Moran 1992; van Buskirk and Steiner 2009; Lande 2014). When environmental variability differs between populations, we can thus expect spatial differences in tolerance to evolve within a species. For example, following the observation that temperature varies more strongly across a day or a season at high than at low latitudes, the climate variability hypothesis attempts to explain latitudinal differences in ectotherms’ thermal tolerance as a consequence of the observed spatial difference in temperature fluctuations (Janzen 1967; Stevens et al. 1989; Gaston and Chown 1999), with higher thermal tolerance selected in more poleward habitats. However, while some empirical studies seem to support the climate variability hypothesis (Addo-Bediako et al. 2000; Vázquez and Stevens 2004), a new study suggests that alternative mechanisms could cause the same spatial differentiation of environmental tolerance, especially range dynamics (Lancaster 2016).

Expansion of a species’ range through adaptation to novel environmental conditions can cause a temporary and local increase in environmental tolerance and phenotypic plasticity, and it potentially creates clinal patterns across species’ ranges. For example, in a recent meta-analysis Lancaster (2016) showed that higher thermal tolerances at high latitudinal margins were found only for insect species that are currently or were recently expanding or shifting their range toward the poles. Instead, ectotherms with stable distributions—mostly endemic or insular species—were shown to have constant thermal tolerance breadths across latitudes. Lancaster (2016) concluded that temporary evolutionary dynamics in the course of range shifts or range expansions are responsible for observed latitudinal clines in thermal tolerance breadth. Such temporary dynamics have been found in analytical models where a transient increase in adaptive phenotypic plasticity appears in populations facing temporal changes in their local environment and thus increase their environmental tolerance to better cope with novel environments (Gavrilets and Scheiner 1993; Lande 2009; Chevin and Lande 2011; Gallet et al. 2014). Following stabilization of the environment, genetic assimilation can cause a reduction in tolerance and phenotypic plasticity and lead to canalization of the genotypes (Waddington 1953; Crispo 2007; Lande 2009; Ergon and Ergon 2016). Similarly, expanding a species’ range by adapting to novel environmental conditions outside the current ecological niche—for instance, by invasive species—can be achieved by a transient increase in plasticity and tolerance in the newly founded populations (Lande 2009, 2015; Chevin and Lande 2011). Plasticity clines can thus have two different origins: range dynamics or clines of climate variability. No attempt has yet been made to distinguish between these two causes and to investigate the effect of range dynamics on tolerance and plasticity evolution in detail.

Species’ range dynamics can result from the colonization of novel habitats by the evolution of a species’ ecological niche and from environmental change. Range expansion through niche expansion is achieved when a species evolves its environmental optimum, its tolerance breadth, or both in novel habitats (Wilson 1961; Thomas et al. 2001; Wiens and Donoghue 2004; Early and Sax 2014; Atwater et al. 2018). Niche evolution is an important driver of invasive species’ range dynamics when the evolution of increased tolerance allows alien species to become invasive in the novel habitat (Brock et al. 2005; Richards et al. 2006; Lande 2009, 2015; Alexander and Edwards 2010; Chevin and Lande 2011). However, the evolution of the ecological niche is not a necessary prerequisite for changing species’ ranges. Large-scale environmental changes can force a species to track its suitable environmental conditions in space and shift its range accordingly. This scenario of range shift may be valid for latitudinal shifts of species’ ranges after the last ice age (Hewitt 1999, 2000) and for more recent responses to global climate warming (Parmesan 2006; Tingley et al. 2009; Talluto et al. 2017). Alternatively, environmental change can act as driver of range expansion when the rate of environmental change is not constant across space, as when global temperature changes much faster at high than at low latitudes (known as polar amplification; see box 5.1 in Stocker et al. 2013). Populations at the poleward range margin may thus follow rapidly shifting local conditions and expand into new geographical areas, whereas trailing-edge populations at the equatorial range margin may face slower local changes and stay in place while keeping the species’ niche constant. Spatial plasticity clines may here differ from scenarios with uniform rates of environmental change or from scenarios of niche evolution in a constant environment. In short, changes in the species’ distribution can be triggered by the evolution of the species’ ecological niche and by environmental change. While plasticity and tolerance evolution as drivers of range expansions have been studied before in two-patch models (Lande 2015), theoretical work that approaches plasticity and tolerance evolution during large-scale environmental change and comprises an entire species’ distribution is missing. Here, a better understanding might not only allow past processes underlying large-scale biogeographic patterns to be unraveled but also allow future evolutionary responses to environmental change to be better predicted.

In this article, we make a distinction between range shift, range expansion, and niche expansion and show that these three evolutionary scenarios differ in the resulting spatial clines of environmental tolerance and phenotypic plasticity. We used two different approaches of modeling environmental tolerance, working with an evolving tolerance curve and evolving norm of reaction (NoR). Using individual-based simulations, we modeled species’ ranges that were limited by external ecological factors, like competitors, diseases, or steep transitions in the environment. We then initiated range dynamics either by shifting these external constraints together with environmental change or by simply withdrawing limitations and allowing for niche expansion. We show that varying histories of environmental change in local populations set on an environmental gradient can act as a driver of tolerance differentiation between populations, even in absence of a species’ niche expansion and spatial differences in environmental variability. We further found that plasticity clines can be in opposite directions depending on whether a species expands its niche into new habitats or follows it across space.

Methods

To model the evolution of environmental tolerance, we used two common approaches: a tolerance curve and an NoR. The tolerance curve describes fitness depending on the environment in a very general way (Lynch and Gabriel 1987), without an environment-dependent phenotypic representation (Lande 2014). Alternatively, modeling phenotypic plasticity using a genotypic reaction norm explicitly maps phenotypes to environments (Via et al. 1995; Whitlock 1996; Chevin et al. 2010; Lande 2014; Svanbäck and Schluter 2012; Valladares et al. 2014). We used individual-based simulations with a modified version of Nemo (Guillaume and Rougemont 2006) to model evolving tolerance curves and reaction norms. The Nemo source code, Nemo init files, and the summary of simulation results are available from the Dryad Digital Repository (https://dx.doi.org/10.5061/dryad.2nc0bn1; Schmid et al. 2019). Individual-based simulations are a powerful tool to study evolution during range dynamics by accounting for genetic drift and evolving additive genetic variances (e.g., Polechová and Barton 2015; Gilbert et al. 2017). Although genetic drift can cause range limits under certain conditions (see Polechová and Barton 2015; Polechová 2018), we chose population sizes large enough to avoid this effect when drift was present but weak.

Life Cycle

In all simulations, dioecious individuals mated within patches at random and without selfing. The fecundity of each female was picked from a Poisson distribution with NoffspringPois(λ) and an average offspring number of λ=5. For each of the females’ offspring, one male was randomly picked from the pool of potential sires, with replacement. Generations were nonoverlapping with even sex ratio on average. Each generation started with breeding (when adults produced offspring), followed by migration (of the offspring), phenotype expression or tolerance curve determination (of the offspring), computation of each offspring’s individual survival probability (depending on the respective environment), selection (stochastic removal of individual offspring depending on their survival probabilities), regulation (all adults died; offspring were discarded randomly until carrying capacity was reached), and aging (offspring was transferred into adult life stage). As a result of this life cycle, plasticity corresponded to developmental or one-shot plasticity (Lande 2015) when phenotypes were expressed once after migration (e.g., seed dispersal) based on the environment of selection with a perfect reliability of the environmental cue, for simplicity. We modeled the connectivity between patches as a stepping stone model with dispersal only between neighboring patches. The species’ range was modeled with absorbing boundaries such that those individuals migrating beyond the (external) niche limit died. Simulations were run with five different migration probabilities (m=0.001, 0.01, 0.1, 0.2, 0.4) and thus cover situations from very low to very high patch connectivity.

Tolerance Curve

We modeled the tolerance curve as a Gaussian function describing fitness (i.e., survival probability) of a genotype depending on the environment (see fig. 1a). The tolerance optimum (t0) and tolerance breadth (t1) were determined by two evolving quantitative traits. We implemented a generalist-specialist trade-off by imposing a constraint between the height and the breadth of the tolerance curve where the maximum fitness Wmax(t1) is a decreasing function of the tolerance breadth (t1; see fig. 1a and eq. [5]). The fitness of an individual with tolerance traits t0 and t1 has fitness W(e) when responding to the environment e:

(1)W(e)=Wmax(t1)exp[(et0)22t1].
Figure 1.
Figure 1.

Two tolerance curves are illustrated in a with identical environmental optima (t0=30) but different tolerance breadths (t1). Given a specialist-generalist trade-off (β>0; see eq. [4]), a higher tolerance translates into a lower maximal fitness. How phenotypic plasticity translates into environmental tolerance is illustrated in b. The two solid lines (light and dark gray) represent two genotypes as linear reaction norms that describe phenotype expression depending on the environment. The two genotypes exhibit different degrees of plasticity (light gray = no plasticity; dark gray = adaptive plasticity). The gray dotted lines show the fitness landscape, with maximum fitness achieved at the black dashed line representing the position of the phenotypic optimum Θ (here assuming Θ=e). The two reaction norms translate into two different tolerance curves depending on the amount of plasticity (dashed lines at the top of the graph). The black dashed line to the right represents the fitness function (i.e., fitness depending on the phenotype) at e=30. The cost of plasticity (β) that reduces the maximal fitness (Wmax) when the absolute value of plasticity deviates from zero is shown in c.

Phenotypic Plasticity

We implemented phenotypic plasticity as a linear NoR (see fig. 1b), as in Schmid and Guillaume (2017; see also Scheiner 1998; Scheiner et al. 2012). The phenotype of each individual (z) is expressed depending on its genotype and the environment in which it develops. The genotype codes for two evolving quantitative traits, the NoR intercept (g0) and the NoR slope (g1). The environment then affects the phenotype depending on the environmental deviation from the reference environment (g2), also called the perception trait (Lande 2009; Ergon and Ergon 2016). We kept the perception trait value constant. The environment-dependent trait value z(e) is then

(2)z(e)=g0+g1(eg2).
The NoR intercept g0 relates to the genotypic value measured in the reference environment e=g2, where the effect of plasticity cancels out (i.e., g1(eg2)=0) and only g0 contributes to the phenotype. The NoR slope g1 controls the degree of plasticity, that is, how strongly the expressed phenotypes differ between environments. Because the environmental position of the perception trait is a key factor in reaction norm evolution (Ergon and Ergon 2016; see also app. A), we ran simulations with three different values of g2 (g2=10, 20, 30).

After the phenotype has been expressed on the basis of the reaction norm, we used a Gaussian selection function to determine the absolute fitness value W(z) as the individual’s survival probability (fig. 1b, black dashed line); W(z) is a Gaussian function of the distance between the expressed phenotype (z(e)) and the phenotypic optimum (Θ), depending on the strength of selection (inversely related to the width of the Gaussian curve ω2) and the cost of plasticity (Wmax(g1)). We assumed that the phenotypic optimum was identical to the environmental value, such that Θ=e and

(3)W(z)=Wmax(g1)exp[(z(e)e)22ω2].

Costs of Plasticity and the Generalist-Specialist Trade-Off

We modeled constitutive costs of plasticity (sensu Chevin et al. 2010) such that Wmax declines with increasing (absolute) values of the NoR slope (g1; fig. 1c) or tolerance breadth (t1). For the cost of plasticity in the NoR model, we used a modified version of the trade-off function from Débarre and Gandon (2010):

(4)Wmax(g1)=(1|αg1|1/β)β,
with the scale parameter α and the shape parameter β, which controls the concavity of the relationship. While α was 0.5 in all of our scenarios, we ran simulations with five different values of β (β=0.0, 0.5, 0.8, 1.0, 1.3). In the absence of costs (β=0.0), maximum fitness is high for all values of plasticity, while β>0 causes a negative relationship between plasticity and maximal fitness, installing a specialist-generalist trade-off. The higher β, the more constrained the evolution of plasticity (fig. 1c). Furthermore, we explored the effect of three different selection strength values (ω2=1, 4, 16) and thus integrated scenarios from very strong to mild stabilizing selection within patches (Kingsolver et al. 2001).

For a better comparison with the evolution of phenotypic plasticity, we translated the cost of plasticity linked to g1 into the cost of tolerance represented by its breadth parameter t1 via g1=1(ω2/t1)1/2 (for the derivation, see app. B), such that

(5)Wmax(t1)=(1|α(1ω2/t1)|1/β)β.

Genetic Parameters

Individuals were diploid in random mating populations. The two evolving quantitative traits were coded on n=20 unlinked quantitative trait loci, each locus contributing additively and pleiotropically to two traits (tolerance curve optimum t0 and breadth t1 or NoR intercept g0 and NoR slope g1). Each allele thus had two additive effects (a0, a1), and the two trait genotypic values were the sum of their respective allelic effects across all loci l and both homolog copies j (e.g., t0=l=120j=12a0lj, t1=l=120j=12a1lj). Mutation effects at pleiotropic loci were drawn from a bivariate normal distribution ηN2(ϵ,M) with average ϵ=(0,0), variances (M1,1,M2,2)0, and zero correlation (M1,2=M2,1=0). Although all loci and mutations were pleiotropic, the absence of mutational correlation ensured that the trait’s genetic covariance remained close to zero on average, allowing the two traits to evolve independently (Chebib and Guillaume 2017). We used a mutation rate of μ=0.0001 per allele and a continuum-of-allele mutation model where mutational effects were added to the existing allelic effects. Mutational variances for phenotypic plasticity were set to M1,1=0.1 for g0 and M2,2=0.001 for g1 in the NoR simulations, such that the phenotypic effects of mutations had a variance of Vm,i=2nμσP2, in the range of Vm0.001, given σP2=M1,1+g12σ2+(eg2)2M2,2 (see app. A). Mutational variances for tolerance curve evolution were set to M1,1=0.1 for t0 and M2,2=1 for t1. We used a higher mutational variance for t1 than for g1, as mutational phenotypic effects σP2 in g1 were environment dependent and increased with the distance between e and g2 along the gradient (see fig. A1). We also ran additional simulations without pleiotropy in which each of the two traits was controlled by 20 separate nonpleiotropic loci. The results obtained were qualitatively very similar to those of the full pleiotropy model (see fig. S1). We thus present only the pleiotropic case in the main text.

Scenarios of Species’ Range Dynamics

Initial Environmental Conditions (Burn-In Simulations)

We modeled 42 habitat patches that were linearly arranged along an environmental gradient connected by nearest-neighbor dispersal (i.e., the stepping stone model). Average environmental values per patch (e¯) ranged from 10 at the left margin to 51 at the right margin with a constant between-patch environmental distance of 1. In the burn-in, we set the initial range within the first 21 patches on the gradient, with average environmental values between 10e¯30 and a carrying capacity of 200 individuals (see fig. 2a, dashed gray line). Patches to the right of the initial range were designated as not habitable and attributed a carrying capacity of zero (e.g., species reached a geographical barrier or the range of a competitor). Environmental conditions also varied randomly within patches across generations and individuals such that the environmental value experienced by each individual (and thus the phenotypic optimum) was picked from a Gaussian distribution with mean e¯ and variance σ2=1 as eN(e¯,1). Consequently, a population experienced within-patch environmental variation resulting either from spatial heterogeneity within a patch or from temporal fluctuations when individuals were born and experienced selection at slightly different points in time. For each parameter combination (table 1) we ran burn-in simulations for 100,000 generations and 20 replicates on a constant average environmental gradient to reach migration-selection-drift balance.

Figure 2.
Figure 2.

Graphs showing the three scenarios of species’ range shifts and environmental change. In a, 42 patches (X-axis) are linearly arranged along an environmental gradient (Y-axis) with environmental values e¯. The dashed line shows the environmental values set in the burn-in simulations before range expansion (21 patches on the left). The solid lines show the environmental values in patches at the end of the expansion/shift simulations after 210 generations for each scenario (range shift, range expansion, and niche expansion; see text for details). Within a single patch, b illustrates the average norm of reaction (NoR; black line), the phenotypic optima (gray dashed line), and the experienced environment before (solid dot) and after (circle) environmental change in the range shift and range expansion scenarios. The lower phenotypic optimum is partially reached thanks to an adaptive plastic response (positive NoR slope). NE = niche expansion; RE = range expansion; RS = range shift.

Table 1.

Explored parameter space of the burn-in and range dynamic simulations (range shift, range expansion, and niche expansion)

ParameterValues
Costs of plasticity (β).0, .5, .8, 1.0, 1.3
Migration rate (m).001, .01, .1, .2, .4
Selection strength (ω2)1, 4, 16
Perception trait value (g2)10, 20, 30

Note. We ran 20 replicates for every parameter combination.

View Table Image

According to Polechová and Barton (2015), we expected no limit to the species’ ranges to evolve in the burn-in populations because genetic drift was not strong enough to impede adaptation at the range margins given our choice of parameters. With an evolving additive genetic variance influenced by gene flow, sharp limits to species’ ranges can evolve from genetic drift only if the environmental gradient (B) is steep and the efficacy of selection relative to the effect of drift from low population sizes (Nσ(s)1/2) is weak, such that B0.15Nσ(s)1/2 (Polechová and Barton 2015). We thus manipulated the range margin conditions to model range shift and expansion scenarios in the follow-up simulations.

Range Shift

After burn-in, the average environmental conditions within patches started to change with rate Δe¯=−0.1 per generation (fig. 2a). As we did not allow the realized ecological niche to evolve in this scenario, patches with environmental values outside the initial species’ range were not available for colonization and their carrying capacities were set to zero. Therefore, patches at the trailing edge became inhabitable when their local environmental value dropped below 9.5. Alternatively, new patches at the front edge became available for colonization when their environmental value reached 30.5. With rate Δe¯=−0.1, the species’ niche shifted completely to the right of the environmental gradient and settled into patches 22–42 in 210 generations.

This scenario applies to species that track their ecological niche in space during environmental change with range margins set by an external factor—for instance, a competitor—and not an intrinsic physiological or evolutionary limit. The species’ realized niche is thus decoupled from its environmental tolerance, and the evolution of environmental tolerance and plasticity has little effect on the position of the ecological margin. Such cases can be detected when a species can establish outside its realized niche in absence of a competitor (Bridle and Vines 2007). We kept the realized ecological niche constant during environmental change for simplicity, such that determinants of external range limit shifted together with the environmental conditions. With this assumption we intended to avoid conflating effects of niche evolution with effects of range shift.

Range Expansion

In this scenario, we explored the consequences of range expansion while maintaining the ecological niche constant. We achieved this by setting variable rates of environmental change across the range, starting with a constant habitable environment at the left margin (constant edge) and a maximum rate of change at the expanding edge of the range (set at Δe¯=−0.1), until patch 42 reached an environmental value of 30.5 and became habitable (fig. 2a). The rate of change in the rest of the range was linearly decreased to maintain a linear environmental gradient among patches as the range increased. In this model, new patches became habitable every 10 generations at the right margin. The time course was similar to the range shift scenario with Δe¯=−0.1, such that patch 42 became habitable after 210 generations of range expansion.

This scenario mimics environmental change with an extreme case of polar amplification (Stocker et al. 2013, box 5.1), when environmental change (e.g., global warming) is stronger at one edge of the gradient (at high latitudes or altitudes) compared with the other edge (at the equator and low altitudes). Again, we assumed the species’ environmental tolerance to be decoupled from external factors determining the ecological margins, and we kept the realized ecological niche constant over time.

Niche Expansion

Finally, we modeled niche expansion on a constant environmental gradient by allowing for range expansion into novel habitats. After the burn-in, the individuals were allowed to colonize the 21 new habitat patches on the right (fig. 2a). To allow for a better comparison with the range shift and range expansion scenarios, the patches on the right were opened for colonization in a stepwise fashion every 10 generations, with the last patch becoming habitable after 210 generations, ensuring that the total time for colonization was the same as with Δe¯=−0.1.

In this scenario, the species was allowed to adapt to novel environmental conditions and to colonize new patches, expanding its range and ecological niche (i.e., to patches with e¯>30). Such a scenario relates, among others, to invasive species that adapt to novel environmental conditions and expand their range.

In all three scenarios, we have run the simulations for an additional 180 generations with stable environmental conditions to study genetic assimilation.

Results

Static Range (Burn-In)

After 100,000 generations in a constant environment, average environmental tolerance (t1) and phenotypic plasticity (g1) evolved to uniform values along the species’ range for most parameter combinations (e.g., see figs. 35, dashed lines). Lower levels of plasticity and tolerance were observed at the range edges compared with the center when migration was high, because genotypes in the edge populations experienced lower environmental variation across generations due to the absorbing boundaries (fig. S2).

Figure 3.
Figure 3.

The norm of reaction (NoR) intercept, the NoR slope, and the resulting phenotype are shown for scenarios with ω2=4, σ2=1, m=0.01, and g2=30 for two different costs (gray: β=0.5; black: β=1.0) before (dashed lines) and after (solid lines) environmental change (range shift [RS], range evolution [RE]; Δe¯=−0.1) or niche evolution (NE).

The average plasticity and tolerance increased with higher migration rate (m), stronger selection (low ω2 values), and lower costs (low β values; fig. 5; table S1). In the absence of costs (β=0), plasticity evolved to be “perfect” (g1=1), and tolerance reached high values (t1>100; fig. S3). The lowest values of plasticity and environmental tolerance in the absence of costs were achieved with the lowest migration rates. The perception trait value (g2) had no effect on the evolved level of plasticity at equilibrium in most of our simulations (fig. S4).

As expected with increasing levels of average plasticity, the genotypic clines in g0 across the range were shallower than the phenotypic clines, which matched the environmental values (figs. 3, S5a). Phenotypic clines across the range were steeper with high phenotypic plasticity (fig. S6). A positive covariance between NoR slope (g1) and NoR intercept (g0) within populations evolved in all three scenarios after 100,000 generations of burn-in, with covariances decreasing toward the range edges (fig. S7).

Environmental Change

We observed the evolution of a spatial cline in tolerance and plasticity in all three scenarios of dynamic species’ ranges, although the cline orientation strongly differed between scenarios. In the range shift scenario, a negative cline evolved with highest plasticity and tolerance at the trailing edge of the distribution (figs. 3a, 4a, S8a, S8d). In the range expansion scenario, plasticity and tolerance were maximized in the middle of the distribution range (figs. 3b, 4b, S8b, S8e), while in the niche expansion scenario a positive cline evolved with the highest values at the expansion front (figs. 3c, 4c, S8c, S8f).

Figure 4.
Figure 4.

The tolerance breadth and the tolerance optimum for scenarios with ω2=4 and σ2=1 for two different costs (gray: β=0.5; black: β=1.0) before (dashed lines) and after (solid lines) environmental change (range shift [RS], range expansion [RE]; Δe¯=−0.1) or niche evolution (NE).

Environmental tolerance and phenotypic plasticity simulations resulted in qualitatively similar patterns, while slight quantitative deviations were observed, when the evolved plasticity clines were steeper than the tolerance clines (compare fig. 3a–3c with fig. 4a–4c). The evolution of steeper clines in plasticity resulted from mutational effects and additive genetic variance that were environment dependent for g1 (and thus varied along the species’ range) but not for t1. Plasticity thus evolved more than tolerance (t1) in environments very different from g2, despite the higher mutational variance in t1.

In line with the expectation of genetic assimilation, the gradient in average tolerance and plasticity leveled out again after range shifts and range expansions (fig. S9), a process that is slower than the rise of spatial differences in t1 and g1 during the range shift (Lande 2009).

The clines in g1 and t1 along the species’ range were steepest, when (1) the strength of selection was high (low ω2; e.g., fig. 5a–5c), (2) costs were high and initial g1 and t1 were small (fig. 5g–5i), and (3) migration rate was low (fig. 5d, 5e). Lowering the migration rate also strongly reduced the rate of range expansion, especially in the niche expansion scenario (fig. 5f). No clines in plasticity or tolerance evolved in the range expansion scenario except for small migration rates (m0.01). Steeper plasticity clines under low migration rates can be explained by a reduced migration load on plasticity evolution and a reduced translocation of genotypes from neighboring populations with adaptive g0 and t0 values. Plasticity clines were shallower in the range expansion scenario than in the range shift scenario because of the smaller extent of environmental change across the range and within patches (see fig. 2; compare fig. 5a, 5d, and 5g with fig. 5b, 5e, and 5h; fig. S10). In comparison, the steepest clines were obtained in the niche expansion scenario (figs. 5c, 5f, 5i, S10), where the species had to adapt to novel environmental conditions for which no adapted genotypes were available within the range.

Figure 5.
Figure 5.

Effects of the strength of selection (ω2; ac), dispersal rate (m; df), and costs (β; gi) on the plasticity clines (g1). Average phenotypic plasticity per patch is given before (dashed lines) and after (solid lines) the range shift scenarios after 210 generations. Unless specified, results are given for high costs (β=1.0), intermediate environmental variance (σ2=1), moderate migration rate (m=0.1) with a perception trait value of g2=30, and strong selection (ω2=1). NE = niche expansion; RE = range expansion; RS = range shift.

In the NoR model of phenotypic plasticity, the position of the perception trait (g2, sometimes also referred to as reference environment) had a strong influence on the spatial variation of the NoR slope (g1). The evolution of tolerance breadth was not affected by this parameter and evolved similar clines as the NoR model when g2=30. Otherwise, plasticity clines were reversed for lower values of g2, revealing the evolution of negative NoR slopes in the range expansion and range shift scenarios (fig. 6a, 6b) but not in the niche expansion scenario (fig. 6c). This maladaptive plasticity did not, however, hinder adaptation to the local conditions and also allowed for phenotypic clines well aligned with the environmental gradient (fig. 6d, 6e). It even favored colonization of new habitats in the niche expansion scenario because moving the reference environment farther left on the range increased the phenotypic effects of allelic variation in plasticity (g1; figs. S11, A1).

Figure 6.
Figure 6.

Phenotypic plasticity (g1; ac) and the corresponding phenotypic values (z; df) are illustrated along the environmental gradient depending on the perception trait value (g2). Simulation results are given for high costs (β=1.0), intermediate environmental variance (σ2=1), low migration rate (m=0.01), and intermediate selection strength (ω2=4). Simulations were run with three distinct g2 values (black: g2=10; dark gray: g2=20; light gray: g2=30). NE = niche expansion; RE = range expansion; RS = range shift.

Discussion

Clines in environmental tolerance and phenotypic plasticity across a species’ range have been reported repeatedly and are considered a critical factor for the persistence of species facing environmental change (Valladares et al. 2014; Bennett et al. 2015). It has been suggested that tolerance clines evolve when species adapt to spatial clines of environmental variability, often resulting in higher tolerance at higher latitude (Addo-Bediako et al. 2000; Sunday et al. 2011) or when species modify their home range during range expansion or range shift (Lancaster 2016). We studied the evolution of tolerance and plasticity clines under the second less well covered hypothesis of species’ range evolution in three scenarios of range shift or expansion that were driven by environmental change or by niche evolution. In the environmental change scenarios, we kept the ecological niche unaltered during environmental change by shifting the external determinants of range limits in parallel with the environment. We showed that in scenarios of range shift or range expansion, without niche evolution higher plasticity and tolerance evolved in the parts of the range that experienced the longest history of environmental change, while lower plasticity was retained in the areas reachable by preadapted genotypes. Therefore, we found the largest plasticity in trailing-edge populations in the range shift scenario and in central populations in the range expansion scenario. In the range shift scenario trailing-edge populations were occupied for the longest time and experienced the largest shift of their local conditions, while in the range expansion scenario the central populations had the most favorable combination of extent and duration of past environmental change, provided migration was sufficiently limited. In contrast, in the scenario of niche evolution during invasion (or range expansion into new environments), the leading-edge populations had the highest plasticity and tolerance because they were colonized by genotypes having to repeatedly adapt to novel environmental conditions. Migration favored plasticity and tolerance during colonization and thus reinforced the evolution of plasticity clines in the niche expansion scenario, while it had an antagonistic effect on cline formation in the range shift and range expansion scenarios because it favored the translocation of preadapted genotypes within the range, limiting the need for a plastic response. In sum, range dynamics can have a profound effect on the evolution of plasticity and tolerance and create clines that resemble empirical patterns. Our results thus confirm that range expansion driven by colonization of novel environments (niche expansion) allows for the evolution of plasticity clines with a spatial increment similar to the pattern found for ectotherms’ latitudinal clines of thermal tolerance (Lancaster 2016), while the other two evolutionary scenarios show results not usually reported in empirical studies. Of course, this does not discard the climate variability hypothesis but shows that an alternative explanation to clinal or latitudinal variation in plasticity and tolerance exists.

Mechanisms Behind the Clines

During all of our range shift and range expansion simulations, plasticity and tolerance evolved in response to novel environmental conditions, experienced over time and across space. Plasticity clines then resulted when genotypes experienced different variability of local environments along the species’ range. However, in contrast to the climate variability hypothesis, the difference in variability was a consequence of the rate of environmental change in the different scenarios and of variation in migration rate instead of seasonal or between-year environmental fluctuations. Interestingly, range dynamics produced plasticity clines similar to those expected under climate variability only in the niche expansion scenario, where the environment was static but genotypes moved along the environmental gradient colonizing new habitats. In that case, plasticity becomes more advantageous at the colonization front, where more plastic or tolerant individuals are selected. Costly plasticity may, however, strongly reduce the pace of colonization.

The evolved clines in environmental tolerance and phenotypic plasticity were transient and leveled out after the stabilization of the species’ distribution (fig. S9), a process known as genetic assimilation (Waddington 1953; Crispo 2007; Lande 2009). However, the decline in plasticity by genetic assimilation happens much slower than the buildup of plasticity (Lande 2009; Scheiner et al. 2017), such that the evolved clines outlasted the actual duration of the range dynamics by far. Our results therefore suggest long-lasting effects of species’ range dynamics on the tolerance and plasticity levels of species, despite their temporary nature.

The movement of genotypes—and thus the rate at which they are exposed to novel environmental conditions across space—was a key component of the evolution of plasticity clines in our simulations. Typically, plasticity was higher with increased migration rates in static environments, which also depended on the steepness of the fitness cline along the environmental gradient. Therefore, the equilibrium plasticity or tolerance level in static environments was a function of the spatial variation of fitness a migrating genotype was exposed to. Strong within-patch stabilizing selection (i.e., strong divergent selection) and high migration rates selected for high equilibrium levels of plasticity (see also Via and Lande 1985; Scheiner 1998; Sultan and Spencer 2002). However, during range evolution migration limited the formation of plasticity clines for two main reasons. First, high migration allowed for the translocation of previously locally adapted genotypes to the corresponding suitable habitat farther away on the shifted environmental gradient, reducing selection pressure on higher plasticity or tolerance. This was particularly the case in the range shift and range expansion scenarios, where plasticity clines were steeper for the lowest migration rates and some plasticity evolved even at the front edge when genotypes could not catch up with their environment. Second, higher migration also imposed an evolutionary load on plasticity by bringing low-plasticity genotypes in environmentally variable populations, which resulted, for instance, in shallower plasticity clines at higher migration in the niche expansion scenario (see fig. 5f).

Empirical Patterns

Our finding of maximized tolerance and plasticity at the trailing edge in the range shift scenario has not, to the best of our knowledge, been described empirically yet. Studies of within-species differences in plasticity or tolerance are rare in general (Valladares et al. 2014), and rear-edge populations of dynamic species ranges (often at the warm margin) are understudied compared with those at the leading edge (Hampe and Petit 2005; Thuiller et al. 2008). This rarity of empirical data is especially unfortunate given current ongoing climate change. While a global temperature rise may progressively favor populations at the cold margin and allow them to expand poleward or upward, tracking favorable conditions (Parmesan 2006; Steinbauer et al. 2018), temperature increases are expected to negatively affect populations at the warm margins because they will experience novel environmental conditions outside the species’ niche (Hampe and Petit 2005; Kremer et al. 2012; Allendorf et al. 2013, pp. 450–451). Thus, warm-edge populations are expected to be under stronger pressure to evolve or plastically respond to climate change (Duputié et al. 2015). Although we have not modeled niche expansion at the trailing edge, we expect evolving phenotypic plasticity to favor niche expansion also at the rear edge of species’ ranges and potentially rescue those populations from extinction.

The only study detecting patterns similar to those derived in our range expansion scenario is, as far as we know, Mägi et al. (2011), who found high morphological plasticity in the distribution center of Agrimonia eupatoria, while the closely related species Agrimonia pilosa at its distribution edge exhibited reduced plasticity for the same traits. However, Mägi et al. (2011) hypothesized that plasticity costs increased with environmental stress level, causing lower plasticity in extreme habitats at the range edges. Our findings in the range shift scenario add another potential explanation for these findings in the context of range expansion under environmental change, and further empirical studies are necessary to discriminate between these alternative hypotheses. The scarcity of empirical evidence for high range-center tolerance and plasticity is not surprising given the evolution of only shallow tolerance and plasticity clines in our simulations (figs. 3, 4).

In contrast to the other two scenarios, our finding of elevated plasticity and tolerance at the leading edge in niche expansion scenario is in line with several other theoretical (Roughgarden 1972; Chevin and Lande 2011; Lande 2015) and empirical (Thomas et al. 2001; Matesanz et al. 2012; Lancaster 2016) studies. Invasive species have been repeatedly found to have expanded their niche in novel habitats by evolving higher phenotypic plasticity and environmental tolerance (Molina-Montenegro and Naya 2012; Atwater et al. 2018). Recently, Lancaster (2016) argued that this process could also explain latitudinal patterns of thermal tolerance in range-expanding ectotherms. Most of the range-expanding species in Lancaster (2016) were invasive species (16 of 20) that expanded from low to high latitudes rather than vice versa. In line with this assumed expansion process, ectotherms were found to “overfill” their cold limit, that is, were found beyond their previously measured cold margin, in a similar meta-analysis (Sunday et al. 2012). Interestingly, within-species increases in thermal tolerance and niche breadth with latitude were observed only at higher latitudes and not (or only weaker) at lower latitudes (Lancaster 2016; Papacostas and Freestone 2016). In our simulations we found constant phenotypic plasticity and environmental tolerance in the part of the species’ range that served as a source for the colonization process, giving further support to the argument of Lancaster (2016).

Comparing the Tolerance Curve and Reaction Norm Approaches

We used two distinct approaches to simulate environmental tolerance evolution during species’ range dynamics—the tolerance curve and the reaction norm—and we observed that they do not always lead to the same qualitative results. Both approaches include two evolving quantitative traits that control the position of the environmental optimum (t0,g0) and the tolerance breadth (t1,g1). Deviations between the two approaches for the same parameter combination arose when the tolerance curve evolved broader tolerance breadths in response to environmental novelty, while the reaction norm approach resulted in maladaptive plasticity (negative g1) and thus smaller tolerance breadths (compare fig. 4a with fig. 6a). Correspondence between tolerance curve and reaction norm evolution depended on the position of the perception trait value g2. The perception trait controls the phenotypic effects of mutations affecting g1, the NoR slope. Following equation (2), to increase the trait value a mutation in g1 must be positive in a patch with e>g2 but negative when e<g2, which can lead to negative slope values to the left of the position of the reference environment e=g2. A negative g1 then represents maladaptive plasticity of a genotype in the patch with environmental condition to which it is adapted. Negative slopes can, however, still be favored when they allow the expression of novel phenotypes during directional selection under environmental shifts. Consequently, the evolution of maladaptive plasticity allowed the phenotypic optimum to be followed over time but came at the expense of maladaptive responses of genotypes to local random environmental fluctuations.

We are aware that these consequences of the perception trait are entirely derived from geometrical reasoning, and a biological explanation of g2 is not immediately obvious. In fact, little attention has been paid to the evolutionary implications of the perception trait, and many theoretical studies simply assumed g2 to be zero, to simplify the reaction norm equation. Studies of genetic canalization more explicitly refer to it as the reference environment where genetic and phenotypic variances are minimized (De Jong 1990; Gavrilets and Scheiner 1993; Lande 2009; Chevin and Lande 2011; see also app. A). However, no model investigated the consequences of varying the reference environments on an environmental gradient, as we did here. One alternative would be to assume that locally adapted populations are canalized in their local environment and thus have each evolved a different value for the perception trait. Ergon and Ergon (2016) extended Lande’s (2009) model for an evolving perception trait and showed that g2 could initially facilitate evolution toward novel phenotypic optima (when it evolves away from the novel e) and subsequently favor canalization (when g2 evolves toward the novel e). Ergon and Ergon (2016) justified their approach by arguing that the perception trait (g2) could be seen as a quantitative trait controlled by gene regulatory processes close to environmental cue perception, like cue activation thresholds for transduction elements. In contrast, the degree of plasticity (g1) might depend more on processes closer to phenotype expression that affect the sensitivity of gene regulation to changing transduction factors. Independently of its underlying genetic basis, empirical evidence of spatial heterogeneity and genetic variation in the perception trait are direly missing.

Model Assumptions

We expect our results to be robust to alternative model assumptions, at least when deviations do not change the cause of spatial or temporal variation in the factors favoring the evolution of clines in plasticity and tolerance. In other words, changes in some of our assumptions (e.g., nonoverlapping generations, developmental plasticity, mutational effect sizes) may affect the steepness of the evolving clines but not their orientation if they do not affect how selection varies across space. However, spatial differences in population sizes, patch connectivity, strength of selection, or cue reliability could also favor the evolution of clines in plasticity and tolerance and would blur the association we observed between tolerance level and range dynamics (e.g., Scheiner 1998; Chevin et al. 2010). Additionally, environmental change and niche expansion might co-occur simultaneously, for instance, when the ecological determinants of range limits do shift at a different speed during environmental change than the focal species. Then spatial differences in tolerance and plasticity might be more complex and result from more than one mechanism.

We have also shown that different genetic architectures of the traits (nonpleiotropic vs. pleiotropic) did not bring large qualitative differences in the range shift scenario, where high trailing-edge tolerance still evolved. Nonpleiotropic architectures mostly enhanced the additive genetic variance of the traits, improving their response to selection in some cases and leading to a stronger effect of the perception trait at high migration rates (see fig. S1). In general, there were no correlated responses to selection caused by pleiotropy in the absence of mutational correlation, although more subtle effects have been reported in a different context (Leimar et al. 2019). Instead, we expect correlated responses with a mutational correlation between NoR intercept and slope (or tolerance optimum and breadth) such that the evolution of tolerance clines might be reversed or counteracted altogether.

Our findings propose an alternative to the environmental variability hypothesis to explain clines in plasticity and tolerance based on a population history of range dynamics and niche evolution. These two hypotheses may prove difficult to disentangle in the wild unless environmental variability and range dynamics can be measured and documented. Then correlative approaches could demonstrate a stronger association between tolerance levels and range position (leading edge, center, trailing edge) than between tolerance and environmental variability levels. Alternatively, comparisons of tolerance clines between dynamic and stationary species’ ranges can be useful, especially when stationary species do not exhibit tolerance clines (Lancaster 2016). However, it still remains to be clarified when homogeneous tolerance and plasticity levels are associated with stationary species’ distributions and when they are not.

We thank Lesley T. Lancaster, Thorbjørn H. Ergon, Luis-Miguel Chevin, Samuel M. Scheiner, Josh van Buskirk, Rassim Khelifa, Tim Connallon, Daniel I. Bolnick, and two anonymous reviewers for lively discussions and valuable comments on the manuscript. Simulations were run on the University of Zurich Science Cloud. M.S. and F.G. were funded by Swiss National Science Foundation grants PP00P3_144846 and PP00P3_176965 to F.G.

Appendix A. Potential Consequences of the Perception Trait

The perception trait could have considerable effects on the evolvability of phenotypic plasticity, as it controls the phenotypic effect sizes of de novo mutations and variance in g1. The perception trait could also lead to the evolution of maladaptive plasticity while adapting to novel environmental conditions. These effects become obvious from the reaction norm equation:

(A1)z(e)=g0+g1(eg2).
A mutation in the reaction norm slope has no phenotypic effect in the environment e=g2 when g1(eg2)=0. In other words, when only the NoR slope would evolve, there would be an invariant point where the reaction norm is fixed and is not affected by changes in g1 (fig. A1). Instead, mutational effects of g1 increase with the distance between e and g2 (fig. A1b). Thus, the perception trait controls the environment dependence of the effect sizes of mutations as well as the additive genetic variance resulting from plasticity.

Based on the assumption that the perception trait itself could evolve, Ergon and Ergon (2016) showed that the perception trait g2 in response to environmental change first evolves away from the novel environmental value e (to allow for an elevated phenotypic variance and mutational effect sizes) and subsequently evolves toward the novel value e (to minimize the effects of [deleterious] mutations when the population reached the phenotypic optimum).

When the evolution of the perception trait is constrained, negative effects of g2 might also result. Given that the perception is fixed to a value smaller than the average experienced environment (e.g., g2=10 with e=20), a decrease in the environment with the associated decrease in the phenotypic optimum (when e=Θ) would favor the evolution of smaller or more negative slope values. While more negative slope values would allow expression of phenotypes more closer to the novel optimum, they would lead to a maladaptive plastic response in the face of further environmental variation (fig. A1c). We therefore argue that the adaptive value of plasticity could differ between the response to directional selection (adaptation to novel environments) and the effects during environmental fluctuations (e.g., temporal fluctuations).

Figure A1.
Figure A1.

These graphs illustrate how the perception trait controls the phenotypic effects of de novo mutations in g1 (a, b) and how a constrained perception trait could cause the evolution of maladaptive plasticity (c). A genotype’s reaction norm (black; g0=20, g1=0.4) and three potential mutant genotypes (blue) after a mutation in g1 (g1=0.1, 0.3, 0.6) with a fixed perception trait of g2=20 is shown in a. In b, the phenotypic variance resulting from mutations in the norm of reaction (NoR) intercept and NoR slope with a continuum-of-allele model is visualized. Mutational effects are picked from a bivariate normal distribution with a variance of M1,1=0.1 for the intercept (g0) and M2,2=0.001 for plasticity (g1) with zero covariances (M1,2=M2,1=0). The results are shown for the analytical solution σP2=M1,1+g12σ2+(eg2)2M2,2 for the following wild-type genotype: g0=20, g1=0, g2=20. In c it is demonstrated how a constrained perception trait could foster the evolution of maladaptive plasticity. The average reaction norm before environmental change (black line) allows expression of a phenotype very close to the phenotypic optimum (black circle) in the environment e=20. The dashed line represents the phenotypic optima depending on the environment. Environmental change to e=17 shifts the phenotypic optimum to 17 (blue circle), such that genotypes with lower (here negative) slope values are favored (e.g., blue line), when the perception trait is fixed to g2=10.

Appendix B. How Plasticity (g1) Translates into Tolerance (t1)

To translate the costs of phenotypic plasticity into costs of environmental tolerance, it is necessary to derive an equation of how g1 translates into t1. As a first step, we search for conditions when fitness depending on reaction norm parameters (g0,g1,g2) equals fitness depending on tolerance curve characteristics (t0,t1):

(B1)exp((g0+g1(eg2)Θ)22ω2)=exp((t0e)22t1).

Multiplying this term by log() and −1 and simplifying the equation by assuming t0=g0=g2=0 as well as e=Θ leads to

(B2)(e(g11))22ω2=e22t1.

After dividing both sides by e2 and multiplying by 2,

(B3)(g11)2ω2=1t1,
it is possible to solve for t1:
(B4)t1=ω2(g11)2.
When solving for g1, the following equation is derived (assuming ω2/t1>0):
(B5)g1=1±ω2/t1.

Literature Cited

References Cited Only in the Supplemental Material

“Then follows a general ‘popular’ account of the forms of apes, their geographical distribution, dwelling-places, food, motions, social life, language, reproduction, education, rearing of young, diseases, life in confinement, and of the apes figured on the Egyptian temples.” Figured: “Tschego ape, profile.” From the review of Brehm’s Animal Life (The American Naturalist, 1877, 11:557–559).

Associate Editor: Tim Connallon

Editor: Daniel I. Bolnick