PNAS PLUS

Pollen dispersal slows geographical range shift and accelerates ecological niche shift under climate change Robin Aguiléea,1, Gaël Raoulb,c, François Roussetd, and Ophélie Ronced a

Laboratoire Évolution et Diversité Biologique, UMR 5174, Université Toulouse III Paul Sabatier, Centre National de la Recherche Scientifique (CNRS), École Nationale de Formation Agronomique (ENFA), 31062 Toulouse, France; bCentre de Mathématiques Appliquées, École Polytechnique, 91128 Palaiseau, France; cCentre d’Écologie Fonctionnelle et Évolutive, UMR 5175, CNRS, 34293 Montpellier, France; and dInstitut des Sciences de l’Évolution, Université de Montpellier, CNRS, Institut de Recherche pour le Développement (IRD), École Pratique des Hautes Études (EPHE), UMR 5554, 34095 Montpellier, France Edited by Nicholas H. Barton, Institute of Science and Technology, Klosterneuburg, and accepted by Editorial Board Member Douglas Futuyma July 13, 2016 (received for review May 12, 2016)

adaptation

| gene flow | spatial heterogeneity | cline | extinction threshold

E

vidence is accumulating that current climate change has already induced extinctions of species (1, 2), and this phenomenon is predicted to accelerate in the next decades (3). In an environment that is changing both in space and in time, species may escape from extinction by changing their geographical distribution such that they track in space the same climatic conditions, by adapting their phenotype to the new climatic conditions where their currently occur, or by a combination of these strategies (4, 5). Poleward and upward altitudinal spatial range shifts of populations have already been observed in many animals and plants (e.g., refs. 1, 6, and 7). Likewise, there is increasing evidence of phenotypic changes induced by climate change, caused by plastic and/or genetic responses (8). The success or failure of these responses strongly depends on dispersal, which influences both colonization and adaptation (3, 9, 10). For plants, the demographic migration (due to seed dispersal) is partially uncoupled from gene flow (due to both pollen and seed dispersal): pollen and seed dispersal are therefore not expected to have the same consequences on spatial range shifts and on ecological niche shifts (9, 11–15). In this paper, we investigate how the dispersal distance of pollen specifically affects the evolutionary and demographic responses of populations under climate change. The challenges associated with adaptation to climate change, and how theoretical studies have attempted to model the responses to these challenges, can be concretely illustrated with the example of Picea sitchensis (Sitka spruce). This tree species naturally occurs on a limited geographical range along the Pacific coast of North America. Bud set date, a phenotype linked with adaptation to climate, is genetically differentiated throughout the range, linearly decreasing from south to north, which suggests divergent selection on this trait across the range, associated with temperature (16, 17). Likewise, many theoretical studies (18–22) www.pnas.org/cgi/doi/10.1073/pnas.1607612113

model divergent selection associated with climatic gradients within the range by assuming that the optimal value of some phenotypic traits varies linearly through space. In these models, gene flow comes predominantly from central populations because of higher population sizes there than at the margins. The strength of the asymmetry of gene flow depends on the steepness of the gradient in population size relative to dispersal distances. Theoretical studies predict that asymmetrical gene flow across steep environmental gradients swamps local adaptation in marginal populations, resulting in the evolution of phenotypic clines flatter than optimal (but see ref. 21) and in the evolution of stable range limits when maladaptation constrains local density (e.g., refs. 14, 23, and 24). Climate change can be modeled by assuming that, as temperature increases throughout the range, the optimal phenotype cline shifts through space. Models then predict that the population escapes from extinction only when the velocity of climate change is below a critical threshold. The population may persist as a traveling wave: the phenotypic cline and density then shift in space at the same speed as climate, tracking with some lag the location where fitness is maximal (18, 21). Although those models allow for in situ adaptation as climate changes, they predict a spatial density shift, but no shift in the fundamental ecological niche, that is, the range of temperatures for which the species has a positive growth rate and the range of phenotypic variation within the whole population remain the same. The original model of Pease et al. (18) assumed that there was a fixed range of environmental conditions (e.g., temperatures) within which adapted populations could grow. Significance Will plants species survive climate warming? By which mechanisms? Ignoring pollen dispersal, previous models have predicted that species may escape extinction by tracking in space the climatic conditions to which they were previously adapted. Using quantitative genetics models, we demonstrate that plants dispersing their pollen could survive under a faster climate change than expected: pollen dispersal slows the range shift in space but facilitates the evolution of new climatic niche limits, such that the species becomes adapted to warmer conditions. Therefore, current predictions of the future spatial distribution of plants, as well as current assessments of the extinction risk due to climate warming, could be significantly altered for plants, especially those dispersing pollen far compared with seeds. Author contributions: O.R. designed research; R.A. and G.R. performed research; and R.A., F.R., and O.R. wrote the paper. The authors declare no conflict of interest. This article is a PNAS Direct Submission. N.H.B. is a Guest Editor invited by the Editorial Board. 1

To whom correspondence should be addressed. Email: [email protected].

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10. 1073/pnas.1607612113/-/DCSupplemental.

PNAS | Published online September 12, 2016 | E5741–E5748

EVOLUTION

Species may survive climate change by migrating to track favorable climates and/or adapting to different climates. Several quantitative genetics models predict that species escaping extinction will change their geographical distribution while keeping the same ecological niche. We introduce pollen dispersal in these models, which affects gene flow but not directly colonization. We show that plant populations may escape extinction because of both spatial range and ecological niche shifts. Exact analytical formulas predict that increasing pollen dispersal distance slows the expected spatial range shift and accelerates the ecological niche shift. There is an optimal distance of pollen dispersal, which maximizes the sustainable rate of climate change. These conclusions hold in simulations relaxing several strong assumptions of our analytical model. Our results imply that, for plants with long distance of pollen dispersal, models assuming niche conservatism may not accurately predict their future distribution under climate change.

Such hard niche limits may explain why species tracked the favorable environment through space rather than evolving a new niche. Duputié et al. (21), however, relaxed this assumption and still reached the same conclusion. Contrarily to these two studies, Polechová et al. (19) assumed that local population growth declined with increasing local density: they then predicted that species would shift their ranges through space faster than climate warming, while evolving a new niche adapted to colder climate. Models including interspecific competition and priority effects reached the reverse conclusion: species would shift their range slower than climate change and adapt to warmer conditions (22). These models overall suggest that any factor slowing or accelerating range shifts may affect the evolution of the fundamental ecological niche under a changing climate. Most of the theoretical studies consider a single dispersal mode (but see refs. 20 and 25). The effect of pollen dispersal on these dynamics in an environment that shifts in time and space is difficult to predict because of its antagonistic effects on adaptation. Pollen dispersal is expected to constrain local adaptation in heterogeneous environments, especially in marginal sink populations, more than seed dispersal does (12–15, 26). The magnitude of this negative effect of pollen dispersal depends on the relative contribution of pollen and seed dispersal, which is known to considerably vary among plant species (27, 28). Conversely, long-distance pollen dispersal could facilitate adaptation to marginal conditions by increasing the amount of local genetic variation and introducing warm-adapted genotypes at the expanding cold margins (9). Whether the positive effect of pollen dispersal will dominate over its detrimental effect is a subject of debate (9, 29) and a source of uncertainty in predicting the responses of plants to climate change (10). Building on previous quantitative genetic models in continuous space (18, 19, 21, 23, 24, 30–32), we analyze here the effect of the distance of pollen dispersal relative to seed dispersal on the

responses of a population to a changing environment. We assume that pollen does not directly limit plant fecundity but that population dynamics depend on local adaptation. We determine critical thresholds allowing a population to avoid extinction, for the velocity of climate change, for the steepness of the environmental gradient, and for the distance of pollen dispersal. Using numerical analyses, we evaluate the robustness of our conclusions with respect to several strong assumptions of our analytical model. Materials and Methods General Assumptions. We consider self-compatible, hermaphroditic plants with no seed bank. Following previous analytical models (18, 19, 21, 23, 24, 30–32), we consider the continuous-time evolution of a population in a continuous space structured by a phenotypic trait, which determines fitness. Population density at time t and at location x is denoted nðx, tÞ; the mean phenotypic trait is denoted
zðx, tÞ (see Table 1 for a summary of the notation). At time t and location x, there is a phenotypic value θðx, tÞ = bðx − vtÞ, which maximizes fitness: this optimal phenotype varies linearly through space, with steepness b, and shifts through space at a constant rate v, which mimics climate change. Without loss of generality, we assume that v is positive. For example, the optimal phenotype θðx, tÞ could be bud set date for Picea sitchensis, which decreases with increasing latitude in the Northern Hemisphere (i.e., from south to north) with slope b < 0 (later bud set in the south, at warmer temperature) (16, 17). The optimal bud set date would shift through space at rate v > 0 because of climate warming, as illustrated in Fig. 1 (dashed lines). We model seed and pollen dispersal as two diffusion processes with coefficient σ 2s and σ 2p, respectively. These processes account for random, unbiased, density independent dispersal according to a Gaussian kernel. Coefficients σ s and σ p measure the root-mean-square distance between the location where seeds and pollen grains are produced and where they, respectively, germinate and fertilize an ovule. We define the parameter γ as the relative contribution of pollen dispersal to total dispersal: γ = ð1=2Þσ 2p =σ 2, where σ 2 = σ 2s + ð1=2Þσ 2p is the mean squared distance of total gene dispersal via seeds and pollen [the factor 1=2 expresses the fact that each haploid pollen grain carries only onehalf the number of genes contained in a seed (13, 34)]. This parameter γ varies with the distance of pollen dispersal relative to the distance of seed

Table 1. Notation, default numerical values, and range of numerical values explored Variables and parameters

Definition; numerical value or expression

x t nðx, tÞ
zðx, tÞ b v θðx, tÞ σs σp σ

Location Time Population size Mean phenotypic trait Slope of the environmental gradient; −0.0324 d·km−1 Rate at which the environmental gradient shifts in space; from 0 to 2,000 km·gen.−1 (i.e., warming up to 5 °C·gen.−1) Optimal phenotype; bðx − vtÞ Mean distance of seed dispersal; 10 km Mean distance of pollen dispersal; fromq0ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi to 447 km (such that γ varies between 0 and 0.999) Mean distance of total gene dispersal; σ 2s + ð1=2Þσ 2p

γ r0 Vp * Vg * Vs

Relative distance of pollen dispersal; ð1=2Þσ 2p =σ 2; from 0 to 0.999 Maximal growth rate of the population; from 1.5 to 2 gen.−1 Phenotypic variance; from 0.1 to 289 d2 Additive genetic variance; such that h2 = Vg =Vp = 0.8; from 0.08 to 231.2 d2 Inverse of the strength of selection; models with constant genetic variance: 5Vp d2·gen.; model with evolving genetic variance: 12.5vLE d2·gen.; this gives similar values of Vs for all models. Scaling factor for carrying capacity; 100 ind. Genetic variance at linkage equilibrium; 0.01–16 d2 Environmental variance; 0.5vLE d2 (such that we observe h2 ≈ 0.8) Per capita birth rate; 3 gen.−1 Scaled genetic potential for adaptation; Vg =ðR0 Vs Þ pffiffiffiffiffiffiffiffi Scaled slope of the environmental gradient; bσ s =ðR0 2Vs Þ pffiffiffi pffiffiffiffiffiffi Scaled rate at which the environmental gradient shifts in space; v 2=ðσ s R0 Þ

k vLE† Ve † β† A B V

The parameter R0 is the scaled intrinsic growth rate: R0 = r0 − Vp =ð2Vs Þ. Numerical values for b, v, and Vp are inspired by the cline in bud set date in Picea sitchensis (16, 17), assuming that current bud set date corresponds to the optimum. Numerical value for σ s and heritability h2 (and therefore Vg) is inspired by Kremer et al. (9). The value of Vs is inspired by Johnson and Barton (33). The numerical value for k is arbitrary (results do not depend on its value). *Parameters specific to models with constant genetic variance (simplified model and model with local density dependence). † Parameters specific to the model with evolving genetic variance.

E5742 | www.pnas.org/cgi/doi/10.1073/pnas.1607612113

Aguilée et al.

PNAS PLUS

asymmetrical gene flow, due to regions of higher density sending more migrants than regions of lower density. The last term in Eq. 3 is the product of the additive genetic variance Vg with the selection gradient ∂
z
r ðx, tÞ, that is, this is the response of the population to selection.

A

Simplified Model. We first consider a simple model allowing for analytical solutions. In this model, we assume global density dependence, meaning that ~ the density experienced R by an individual at location x, nðx, tÞ, equals total population size λðtÞd nðx, tÞdx. This means that the individual growth rate is depressed by all other individuals independently of their location, and this implies that the decrease of the individual growth rate due to density dependence is the same for all individuals. We will look for equilibria with constant total population size, so we assume that λðtÞ does actually not depend on time. We also assume here that the additive genetic variance Vg is constant in space and in time. To analyze this model, it is helpful to scale variables and parameters, as done in Kirkpatrick and Barton (24) (from now, we also drop the mention that variables depend on time and space to simplify the notation). In particular, we scale time as T = R0 t, where R0 = r0 − Vp =ð2Vs Þ is the maximum growth rate of a population with mean phenotype matching the optimal phenotype. Carrying capacity is scaled by relative growth rate as K = kR0 =r0. Space is scaled by pffiffiffiffiffiffiffiffi the seed dispersal distance and scaled time (X = x 2R0 =σ s), population size by the scaled carrying capacity (N = n=K and Λ = λ=K), and mean phenotype by pffiffiffiffiffiffiffiffiffiffi the strength of selection (Z =
z= R0 Vs). Note that this rescaling uses the dispersal distance of seeds only. This is convenient for analyzing the effect of the distance of pollen dispersal at a fixed distance of seed dispersal. After rescaling, the model formed by Eqs. 1 and 3 can be written with only four parameters:

Fig. 1. Traveling-wave equilibrium (Eq. 6). (A) Scaled population density N (solid black curve), scaled mean phenotype Z (solid black line), and scaled mean growth rate R (solid gray curve; Eq. 5) as a function of space at the travelingwave equilibrium. The dashed line corresponds to the environmental gradient. The thick segment on the y axis depicts the realized ecological niche. (B) State of the population 10 time units later. Parameters of the dynamic equilibrium are depicted (Eqs. 6–11). Arrows have a length proportional to the rate of the shifts they depict. Parameter values: A = 0.1, B = −0.05, V = 2, γ = 0.65.

dispersal: γ = 0 when only seeds disperse and pollen grains fertilize ovules locally; γ becomes close to 1 when pollen dispersal distance becomes large compared with seed dispersal distance. Change in population size is described by the following: ∂t nðx, tÞ =

σ 2s ∂x,x nðx, tÞ + n
r ðx, tÞ, 2

[1]

where ∂x,x denotes the second-order partial derivative with respect to space and
r ðx, tÞ is the intrinsic growth rate averaged over all phenotypes at location x and at time t. The first term of the right-hand side of Eq. 1 describes the effect of seed dispersal. Pollen dispersal does not contribute here to change in population size because we assume that pollen is produced in sufficient amount so that all ovules are fertilized. The second term, detailed below, accounts for the birth and death of individuals. ~ ðx, tÞ=kÞ, The fitness of an individual with an optimal phenotype is r0 ð1 − n where the factor k scales local carrying capacity and where the expression of ~ ðx, tÞ, will vary depending on the local density experienced by an individual, n assumptions (see subsequent sections). Individual fitness rðx, t, zÞ decreases quadratically as the individual phenotype z deviates from the optimum θðx, tÞ: ~ ðx, tÞ=kÞ − ðz − θðx, tÞÞ2 =ð2Vs Þ, with 1=Vs being the strength rðx, t, zÞ = r0 ð1 − n of selection. Integrating rðx, t, zÞ over the phenotypic distribution with variance Vp gives the mean growth rate of a population with mean phenotype
z:   ~ ðx, tÞ n ð
zðx, tÞ − θðx, tÞÞ2 Vp
r ðx, tÞ = r0 1 − − − . k 2Vs 2Vs

[2]

The second term of the right-hand side of Eq. 2 is the lag load (35). The last term is the loss of mean fitness resulting from variation of the phenotypes around the mean. Change in mean phenotype is expressed as follows: ∂t
zðx, tÞ =

σ2 ∂x,x
zðx, tÞ + σ 2 ∂x
zðx, tÞ  ∂x logðnðx, tÞÞ + Vg ∂
z
r ðx, tÞ. 2

[3]

The first term of the right-hand side of Eq. 3 corresponds to the random dispersal of genes carried by pollen and seeds. The second term results from

Aguilée et al.

∂T N = ∂X,X N + NR, ∂T Z =

1 2 ∂X,X Z + ∂X Z  ∂X logðNÞ − AðZ − BðX − VT ÞÞ, 1−γ 1−γ

[4]

where 1 R = 1 − Λ − ðZ − BðX − VTÞÞ2 2

[5]

is the scaled mean growth rate at location X and at time T. In Eq. 4, the parameter A measures the genetic potential for adaptation to the environmental gradient: A = Vg =ðR0 Vs Þ. The parameter B is the slope of the pffiffiffiffiffiffiffiffi scaled environmental gradient: B = bσ s =ðR0 2Vs Þ. Finally, Vpis scaled rate ffiffiffi thep ffiffiffiffiffiffi at which the environmental gradient shifts in space: V = v 2=ðσ s R0 Þ. Note that B and V depend on the dispersal distance of seeds. The fourth parameter is γ, the distance of pollen dispersal relative to the distance of total dispersal (General Assumptions). We analytically identified different equilibria for this model and the conditions for their existence. Stability of the equilibria was investigated by numerically solving the model with various initial conditions (Appendix S1. Numerical Resolution of Eqs. 1 and 3). Model with Local Density Dependence. We numerically solve Eqs. 1 and 3 assuming local density dependence, that is, the density experienced by an in~ ðx, tÞ, equals local population size nðx, tÞ (Appendix S1. dividual at location x, n Numerical Resolution of Eqs. 1 and 3). By contrast to the simplified model, the individual growth rate is here depressed only by individuals at the same location, which implies that the decrease of the individual growth rate may not be the same for all individuals. We use parameter values inspired by the cline in bud set date in Picea sitchensis (16, 17) assuming that current bud set date corresponds to the optimum, and by Johnson and Barton (33) and Kremer et al. (9) for parameters whose value is unknown for bud set date in Picea sitchensis (Table 1). Model with Evolving Genetic Variance. We relaxed the assumption of a genetic variance Vg constant through space and time using a third model describing the evolution of the distribution of genotypes at each location (Appendix S2. Model with Evolving Genetic Variance). To describe the genetic transmission of the trait, we use here the infinitesimal model, where the genetic variance at linkage equilibrium, vLE, remains constant (36). The full distribution of offspring (i.e., in the population, unconditional to the parental genotypes) is not constrained to be Gaussian, and the additive genetic variance in the population can change both in space and in time due to changes in linkage disequilibrium. We further assume that density dependence is local and independent of the genotype. We solve this model numerically using the parameter from Table 1.

PNAS | Published online September 12, 2016 | E5743

EVOLUTION

B

Results For all models, we found that, if the population manages to avoid extinction, it persists at the equilibrium either as a traveling wave tracking the shifting environmental gradient because of both a spatial density shift and an ecological niche shift, or as a population with a uniform density on the whole available space. We first detail the effects of pollen dispersal on the properties and existence of each equilibrium in the simplified model. Then we evaluate the robustness of these effects using the results of the two other models.

A

Traveling-Wave Equilibrium. The population may reach a dynamic

equilibrium where we have the following:

! ðX − Cn T − Ln Þ2 NðX, TÞ = N0 exp − , 2Vn

[6]

ZðX, TÞ = SðX − Cn T − Ln Þ + DT (all statements in this section are demonstrated in Appendix S3. Resolution of the Simplified Model: Traveling-Wave Equilibrium). Density pffiffiffiffiffiffi N is a Gaussian traveling wave with constant relative width Vn. Under a shifting environmental gradient (V > 0), the center of the range, where population size is maximal, lags behind the locality where the mean growth rate R is maximal, that is, where perfectly adapted individuals are located (Fig. 1). This spatial lag is proportional to another spatial lag, the constant distance −Ln that the population would have to move for individuals at the center of the range to reach the place where they would be perfectly adapted (Fig. 1). The mean phenotypic trait Z is a linear cline with slope S. The population density moves across space at a constant rate Cn. The parameter D in Eq. 6 is the rate of phenotypic change of individuals at the center ofpthe ffiffiffiffiffiffi range (i.e., at X = Cn T + Ln). Because the size of the range Vn remains constant through time, the phenotypic range of living individuals, that is, the whole realized ecological niche of the population, shifts at this rate D (Fig. 1), which is therefore a measure of the rate of ecological niche shift. We will show that, with pollen dispersal, the population density shifts in space slower than the environment and the ecological niche shifts in the direction of the environmental change. Hereafter, we highlight the effects of pollen dispersal on the traveling-wave equilibrium. The effects of the other parameters are detailed in Appendix S4. Traveling-Wave Equilibrium: Other Parameters Effect. Appendix S5. Traveling-Wave Equilibrium: Comparison with Duputié et al. (21) shows that all our results converge to those in Duputié et al. (21) [the results of which are themselves consistent with those of Pease et al. (18) and others] when pollen does not disperse. Pollen dispersal flattens the phenotypic cline. At the traveling-wave equilibrium, the slope of the phenotypic cline equals the following: A S = signðBÞ pffiffiffi ð1 − γÞ. 2

[7]

Increasing the pollen dispersal distance (i.e., increasing γ from 0 to 1) thus flattens the phenotypic cline (Fig. 2A). Pollen dispersal indeed induces a genetic swamping effect so that individuals in the whole range of the population become phenotypically more similar to each other when the distance of pollen dispersal increases. Pollen dispersal contracts the spatial range. The size of the spatial range of the population, measured by Vn, is reduced by pollen dispersal (Fig. 2B), as predicted by the following: Vn =

2

pffiffiffi . jBj 2 − Að1 − γÞ

[8]

A flatter phenotypic cline as pollen disperses further (Eq. 7) results necessarily in maladaptation increasing faster with distance away from the location where individuals are perfectly adapted E5744 | www.pnas.org/cgi/doi/10.1073/pnas.1607612113

B

C

Fig. 2. Effect of the distance of pollen dispersal (γ, black vs. gray elements) on the properties of the scaled mean phenotype Z (A), of the scaled population density N (B), and of the scaled mean growth rate R (C) at the traveling-wave equilibrium (Eqs. 5–11). The dashed line on A corresponds to the environmental gradient. Arrows have a length proportional to the rates they depict (see text). Parameter values: A = 0.1, B = −0.05, V = 1.

[Fig. 2A: maladaptation is the mismatch between realized mean phenotype (solid line) and optimal phenotype (dashed line)], which is mirrored by the mean growth rate of the population decreasing faster with distance from the location with maximal fitness (Fig. 2C). Therefore, individuals having a positive growth rate occupy a narrower range (Fig. 2C), which contracts the spatial range of the population. Pollen dispersal slows down the spatial density wave. When the distance of pollen dispersal increases, our model predicts that the rate Cn at which the population density shifts across space decreases: Cn =

V . 1 + jBjApffiffi2 γ

[9]

Indeed, when the phenotypic cline is flatter (i.e., with pollen dispersing further), the location with highest mean growth rate is closer to the center of the range (Fig. 2 B and C). This location determines where most individuals will be in the future, and Aguilée et al.

[10]

This ecological niche shift occurs in the direction of climate change (Fig. 1): in the context of climate warming, the species niche would phenotypically shift to be adapted to warmer temperature and would lose the ability to grow under the previously colder temperature where it occurred. When pollen disperses, the density moves in space slower than the environment. To have avoided extinction, the population must then have adapted to the new conditions that it experiences throughout its range, that is, the ecological niche must shift in the direction of climate change. The further pollen is dispersed, the slower the spatial range shift (Eq. 9), and therefore the greater the ecological niche shift (i.e., jDj increases with γ increasing; Eq. 10; Fig. 2A). When pollen does not disperse, the density moves in space at the same rate as the environment; there is thus no ecological niche shift (D = 0 when γ = 0; Eq. 10). Pollen dispersal reduces the spatial lag of the phenotypic cline. As a result of the flatter phenotypic cline (Eq. 7) and of the greater ecological niche shift in the direction of climate change (Eq. 10), the phenotypic cline tracks the environmental gradient more closely when the distance of pollen dispersal increases (i.e., jLn j decreases as γ increases; Fig. 2A): Ln = −

V pffiffiffi . jBj 2 + Aγ

[11]

Existence and stability of the traveling-wave equilibrium. The population may escape from extinction as a traveling wave under two conditions (Appendix S3. Resolution of the Simplified Model: TravelingWave Equilibrium). First, the rate V of the environmental shift must be below the threshold: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi pffiffiffi   A jBj 2 A crit 1− + ð1 − γÞ. [12] VTW = 2 1 + pffiffiffi γ 2 2 jBj 2 crit If the environment shifts faster than VTW , then the population goes extinct because it fails to shift rapidly enough its spatial distribution and/or its ecological niche. crit Pollen dispersal has antagonistic effects on the critical rate VTW : increasing the distance of pollen dispersal reduces the spatial lag jLn j of the phenotypic cline (Eq. 11) so that the population may persist under a faster environmental shift, but higher pollen dispersal contracts the range of the population (Eq. 8), which compromises its ability to track a faster environmental shift (Fig. 3). As a result of these antagonistic effects, there is an optimal distance of pollen dispersal (relative to the distance of seed dispersal), γ opt TW, which maximizes the critical rate of sustainable change (Fig. 3):   pffiffiffi 4 2 1 opt jBj 2 − . [13] γ TW = − 3 A 3

The second condition for the population to persist as a traveling wave is that the slope of the environmental gradient (in absolute value, jBj) is above the threshold: A jBcrit TW j = pffiffiffi ð1 − γÞ. 2

[14]

The population can reach the traveling-wave equilibrium for shallower environmental gradients when pollen disperses further. Aguilée et al.

B

C

Fig. 3. Scaled maximal sustainable rate of the environmental shift for the crit traveling-wave equilibrium (solid curves, VTW , Eq. 12) and for the uniformcrit , Eq. 16) as a function of the distance of density equilibrium (dashed line, VUD pollen dispersal (γ), for three values of the slope B of the environmental gradient (A–C). The vertical dotted lines indicate the critical distance of pollen dispersal (γ crit TW, Eq. 15). Letters E, T, U, and B indicate which equilibrium exists in each part of the parameter space. We used A = 2. For B = −1 (B), the crit optimal distance of pollen dispersal γ opt TW (Eq. 13) is 0.63 (i.e., VTW is maximal opt for γ = 0.63). For B = −0.4 (A), γ opt TW = 1; for B = −2 (C), γ TW = 0.

Eq. 14 can be written as a critical distance of pollen dispersal γ crit TW above which the population may reach the traveling-wave equilibrium (Fig. 3): γ crit TW

pffiffiffi jBj 2 . =1− A

[15]

Numerical resolution of the simplified model revealed that the traveling-wave equilibrium is locally stable when it exists, except when the uniform-density equilibrium (see next section) also exists and γ is close to the threshold γ crit TW (Fig. 4A). Uniform-Density Equilibrium. The population may reach an alternative equilibrium where the population density is constant across space and the mean phenotype is a linear cline with a slope matching that of the optimal phenotype (Appendix S6. Resolution of the Simplified Model: Uniform-Density Equilibrium). The persistence of the population at this equilibrium is possible only if the rate of environmental change is below the threshold: pffiffiffi A 2 crit . [16] = VUD jBj

This equilibrium was found in previous studies without pollen dispersal (19, 24). The condition for its existence is independent of the distance of pollen dispersal (Eq. 16). Numerical resolution revealed that, when this equilibrium exists, it is always locally stable (Fig. 4A). Robustness of the Results. Results with local density dependence. Now assuming that density

dependence is local (Materials and Methods, Model with Local Density Dependence and Appendix S7. Robustness of the Results, PNAS | Published online September 12, 2016 | E5745

EVOLUTION

ABV γ . D = − pffiffiffi jBj 2 + Aγ

A

PNAS PLUS

therefore sets the speed of the traveling wave. The shift of the density occurs at the same rate as the environment only when pollen does not disperse (Cn = V when γ = 0; Eq. 9); otherwise, it shifts more slowly than the environment (Cn < V when γ > 0; Eq. 9; Fig. 1). Pollen dispersal entails adaptation to changing local condition. At the traveling-wave equilibrium, the ecological niche shifts at a constant rate D:

A

B

C

Fig. 4. Local stability of the equilibria of the three models (A–C), evaluated by numerical resolution (Appendix S1. Numerical Resolution of Eqs. 1 and 3 and Appendix S2. Model with Evolving Genetic Variance). Symbols filling the graphs crit crit , vUD , and γ crit show the observed equilibria. Lines depict the thresholds vTW TW using nonscaled variables (Table S1). For the model with evolving genetic variance, we first determined the predicted possible equilibria with the expressions of nonscaled thresholds (Table S1), using the observed value of the genetic variance at the end of the numerical resolution. Then, we determined the crit crit observed thresholds vTW , vUD , and γ crit TW as the values of v and γ such that we predict the appropriate switch of equilibria given the observed value of the genetic variance. Parameter values: default parameter values (Table 1) and r0 = 2 gen.−1. For the model with evolving genetic variance, vLE = 1 d2. For the two other models, Vp = 2.5 d2 and Vg = 2 d2. Note that with these parameter values, γ crit TW > 0: the traveling-wave equilibrium therefore does not exist at γ = 0.

S7.1. Results with Local Density Dependence), we found that the crit thresholds VTW and γ crit TW are well approximated by Eqs. 12 and crit 15, respectively, whereas Eq. 16 slightly underestimates VUD (Fig. 4B). When both the uniform-density and the traveling-wave equilibria exist, the traveling-wave equilibrium may be reached only for very high values of the distance of pollen dispersal (Fig. 4B). Similarly to the simplified model, at the traveling-wave equilibrium, with local density dependence, increasing the distance of pollen dispersal flattens the phenotypic cline, contracts the spatial range, slows the spatial density shift, reduces the spatial lag of the phenotypic cline, and accelerates the ecological niche shift (Fig. 5). Local density dependence has nevertheless noticeable effects on the properties of the traveling-wave equilibrium (Fig. 5), in particular accelerating the spatial shift in density, widening the range, and modifying the rate or even the direction of niche shift. The effects are, however, weak for parameter values close to the extinction crit thresholds VTW and γ crit TW (Appendix S7. Robustness of the Results, S7.1. Results with Local Density Dependence and Fig. S1). Results with evolving genetic variance. In this model, density dependence is local and the additive genetic variance can change both in space and in time due to changes in linkage disequilibrium (Materials and Methods, Model with Evolving Genetic Variance, and Appendix S7. Robustness of the Results, S7.2. Results with Evolving Genetic Variance and Fig. S2). The existence and local stability of both the traveling-wave and the uniform-density equilibria are rather similar to the results found with the other models (Fig. 4C). Nevertheless, the critical rate of the environcrit mental shift VTW is underestimated at low distance of pollen dispersal and stops varying when the distance of pollen dispersal E5746 | www.pnas.org/cgi/doi/10.1073/pnas.1607612113

becomes very large (Fig. 4C). Despite a slight decrease of the genetic variance as the distance of pollen dispersal increases (Appendix S7. Robustness of the Results, S7.2. Results with Evolving Genetic Variance and Fig. S3), Fig. 5 shows that allowing the genetic variance to evolve does not qualitatively alter the effects of the distance of pollen dispersal on the features of the travelingwave equilibrium, and it has a weak quantitative effect, much weaker than the effect of considering local density dependence vs. global density dependence. Discussion For a population living in a spatially heterogeneous landscape and experiencing climate change, the original model by Pease et al. (18) predicted that a population escaping extinction shifts its spatial range but keeps the same fundamental niche limits (see also ref. 21). Even though the model allowed for genetic evolution, no ecological niche shift was predicted when niche limits arise from the genetic swamping effect of dispersal. Under the same assumptions, we showed here that plants species surviving climate change are expected simultaneously to shift their geographical range and to shift their ecological niche in the

A

B

C

D

E

Fig. 5. Features of the traveling-wave equilibrium (A–E) as a function of pollen dispersal distance (γ), using nonscaled variables (Table S1). On C, the dashed horizontal line indicates the speed of environmental change. On E, the dashed horizontal line indicates no ecological niche shift. Parameter values: default parameter values (Table 1) and r0 = 1.5 gen.−1, v = 11 km·gen.−1. For the model with evolving genetic variance, we used vLE = 0.05 d2. For the two other models, we used Vp = 0.125 d2 and Vg = 0.1 d2, which is the observed value of the additive genetic variance at equilibrium with the model with evolving genetic variance at γ = 0. The difference between the model with local density dependence (circles) and the model with evolving genetic variance (triangles) therefore depends on the effect of the distance of pollen dispersal on the equilibrium value of the additive genetic variance.

Aguilée et al.

1. Parmesan C (2006) Ecological and evolutionary responses to recent climate change. Annu Rev Ecol Evol Syst 37(1):637–669. 2. Selwood KE, McGeoch MA, Mac-Nally R (2015) The effects of climate change and land-use change on demographic rates and population viability. Biol Rev Camb Philos Soc 90(3):837–853. 3. Urban MC (2015) Accelerating extinction risk from climate change. Nature 348(6234):571–573. 4. Davis MB, Shaw RG, Etterson JR (2005) Evolutionary responses to changing climate. Ecology 86(7):1704–1714.

5. Bellard C, Bertelsmeier C, Leadley P, Thuiller W, Courchamp F (2012) Impacts of climate change on the future of biodiversity. Ecol Lett 15(4):365–377. 6. Chen IC, Hill JK, Ohlemüller R, Roy DB, Thomas CD (2011) Rapid range shifts of species associated with high levels of climate warming. Science 333(6045):1024–1026. 7. Hill JK, Griffiths HM, Thomas CD (2011) Climate change and evolutionary adaptations at species’ range margins. Annu Rev Entomol 56(1):143–159. 8. Merilä J, Hendry AP (2014) Climate change, adaptation, and phenotypic plasticity: The problem and the evidence. Evol Appl 7(1):1–14.

Aguilée et al.

ACKNOWLEDGMENTS. We thank A. Duputié, D. Waxman, and S. Mirrahimi for help on the analysis of the models. We are grateful to two anonymous reviewers for comments and suggestions that improved our manuscript. We thank P. Solbès for support on the EDB-Calc system. R.A. was partly supported by the How Does Evolution Affect Extinction And Species Range Dynamics (EVORANGE) (ANR-09-PEXT-01102) project from the French “Agence Nationale de la Recherche” (allocated to O.R.). This work was supported by the Mechanisms of Adaptation to Climate Change (MeCC) (ANR-13-ADAPT-0006) and How Does Selfing Affect Adaptation (SEAD) (ANR-13-ADAPT-0011) projects (allocated to O.R.), by the Modèles Mathématiques pour la Biologie Évolutive (MODEVOL) project (ANR-13-JS01-0009) (allocated to G.R.), and by the French Laboratory of Excellence project “TULIP” (ANR-10-LABX-41). This work was performed using high-performance computing resources from Calcul en Midi-Pyrénées [CALMIP, Grant 2013-153 (allocated to R.A.)], using the cluster of the Institut des Sciences de l’Évolution, and using the cluster of the laboratory Évolution et Diversité Biologique (EDB-Calc, which includes software developed by the Rocks Cluster Group at the San Diego Supercomputer Center, University of California).

PNAS | Published online September 12, 2016 | E5747

PNAS PLUS

dispersal, we found different results, even including the opposite one where the traveling-wave equilibrium replaces the uniform-density equilibrium when the genetic variance evolves (Fig. 4 B and C). Then, the genetic swamping effect of pollen dispersal reducing local adaptation and therefore reducing genetic variance between populations dominates over the input of alleles from distant populations, which increases local genetic variance. For all of the parameter values we tested, pollen dispersal thus helps adapting to a changing climate because it slows the spatial range shift but not because of an increase of local genetic variance. This conclusion holds when traits subject to selection are determined by many independent, additive loci, with evolution of the genetic variance due to changes in linkage disequilibrium [i.e., the infinitesimal model (36)]. Although this model of inheritance approximates well more realistic explicit multilocus models (43, 44), more investigation would be valuable to know how other genetic architectures, in particular those allowing for the evolution of genetic variance at linkage equilibrium, would alter our conclusion. Other mechanisms than pollen dispersal have been described as generating species range limits and causing ecological niche shifts under climate warming. Interspecific competition may prevent successful migration toward colder climate and result in the evolution of the ecological niche, such that species become adapted to warmer temperature when dispersal is limited by habitat fragmentation (22, 45, 46). Atkins and Travis (20) described similar space-blocking effects, within a species, slowing spatial range shifts when competition involves priority effects. Conversely, Polechová et al. (19) showed that local density dependence may accelerate the speed of the traveling wave, resulting in species range moving faster than climatic change and leading to an ecological niche shift in the opposite direction of climate change, that is, such that species adapt at its expanding edge to colder temperatures than temperatures previously included in its niche. We here found a similar effect of local density dependence when pollen dispersal is low (Fig. 5 C and E). Similarly, selection of higher dispersal abilities at range margins may speed up spatial range shifts (22, 47), helping surviving climate change with ecological niche shift in the opposite direction to climate warming. Note that a lack of evolutionary potential at some key ecological traits would prevent any ecological niche shift, even in the presence of one or several of the above mechanisms (48). Our study thus contributes to a set of ecoevolutionary models predicting that the direction and magnitude of ecological niche shifts and the speed of spatial range shifts may be highly variable between species, consistently with empirical findings of various speeds of range shifts associated to current climate warming (6).

EVOLUTION

direction of environmental change when pollen dispersal distance is large enough (Fig. 1). The genetic swamping effect of pollen dispersal restricting local adaptation (12–15, 25) indeed contracts the geographical range of the population and prevents the population from tracking climate in space. Therefore, a population avoiding extinction necessarily adapts locally to the new conditions that it experiences throughout its range, that is, the population shifts its ecological niche in the direction of climate change. We showed that ecological niche shift increasingly dominates spatial range shift when pollen disperses increasingly far compared with seeds (Fig. 2). We demonstrated that the qualitative effect of pollen dispersal is independent of the strongest assumptions of our analytical model (global density dependence and constant genetic variance; Figs. 4 and 5; assumptions on density dependence have nevertheless noticeable consequences, as detailed below). A major consequence of these results is that ecological niche shift cannot be a priori omitted when forecasting future spatial distribution of species. Studies aiming at predicting the future spatial distribution of species use mostly correlative methods that determine the locations where the future climatic conditions will be similar to the conditions where a species currently occur (37). These methods have become more and more refined in recent years, with the addition for example of dispersal limitations and explicit demographic processes (38) or of biotic interactions (39), but they still assume that the fundamental ecological niche remains constant (but see ref. 40). Our results emphasize that considering the chance for adapting to new climatic conditions is actually essential for plants, especially for plants whose pollen disperses far compared with seeds. We expect this conclusion to apply also to animals with separate sexes in which only females (analogous to seeds) directly contribute to population growth, and male (analogous to pollen) dispersal is substantially greater than female dispersal (15, 41). Another crucial consequence of the ecological niche shift predicted by our models is that pollen dispersal may allow plants to survive faster climate change than without pollen dispersal (Figs. 3 and 4, and Eq. 12). We indeed showed (Eq. 13) that, unless the environment is spatially very heterogeneous, the direct maladaptive effect of pollen dispersal (genetic swamping) may be overcome by its indirect beneficial consequence (ecological niche shift). Similarly, Atkins and Travis (20) found that the risk of extinction under climate change is higher with shorter pollen dispersal distance. As a result, the maximal sustainable velocity of climate change may have been previously underestimated for plants, in particular those with a long distance of pollen dispersal compared with that of seeds, as is expected for different tree species for which the maximal distance of pollen dispersal is several orders of magnitude larger than that of seeds (9). Although reaching contrasted conclusions, several reviews of the literature (9, 10, 29) highlighted the possible beneficial effect of pollen dispersal in a context of climate change because it may provide the local genetic diversity necessary to adapt at the leading edge of the population. This possibility is supported in particular by Barton (31) and Polechová et al. (19), who showed with models similar to ours, but without pollen dispersal, that, when the genetic variance is allowed to evolve, the uniform-density equilibrium (i.e., a population with reasonable adaptation on an unlimited geographical range) replaces the traveling-wave equilibrium (i.e., a population with restricted adaptation on a limited geographical range), unless the population is subject to strong genetic drift (14, 42). With pollen

9. Kremer A, et al. (2012) Long distance gene flow and adaptation of forest trees to rapid climate change. Ecol Lett 15(4):378–392. 10. Corlett RT, Westcott DA (2013) Will plant movements keep up with climate change? Trends Ecol Evol 28(8):482–488. 11. Thrall PH, Richards CM, McCauley DE, Antonovics J (1998) Metapopulation collapse: The consequences of limited gene-flow in spatially structured populations. Modeling Spatiotemporal Dynamics in Ecology, eds Bascompte J, Sole RV (Springer, New York). 12. Butlin RK, Bridle JR, Kawata M (2003) Genetics and the boundaries of species’ distributions. Macroecology: Concepts and Consequences, eds Blackburn TM, Gaston KJ (Cambridge Univ Press, Cambridge, UK), pp 274–295. 13. Hu XS, He F (2006) Seed and pollen flow in expanding a species’ range. J Theor Biol 240(4):662–672. 14. Bridle JR, Polechová J, Kawata M, Butlin RK (2010) Why is adaptation prevented at ecological margins? New insights from individual-based simulations. Ecol Lett 13(4): 485–494. 15. Aguilée R, Shaw FH, Rousset F, Shaw RG, Ronce O (2013) How does pollen versus seed dispersal affect niche evolution? Evolution 67(3):792–805. 16. Mimura M, Aitken SN (2007) Adaptive gradients and isolation-by-distance with postglacial migration in Picea sitchensis. Heredity 99(2):224–232. 17. Aitken SN, Yeaman S, Holliday JA, Wang T, Curtis-McLane S (2008) Adaptation, migration or extirpation: Climate change outcomes for tree populations. Evol Appl 1(1): 95–111. 18. Pease CM, Lande R, Bull JJ (1989) A model of population growth, dispersal and evolution in a changing environment. Ecology 70(6):1657–1664. 19. Polechová J, Barton N, Marion G (2009) Species’ range: Adaptation in space and time. Am Nat 174(5):E186–E204. 20. Atkins KE, Travis JMJ (2010) Local adaptation and the evolution of species’ ranges under climate change. J Theor Biol 266(3):449–457. 21. Duputié A, Massol F, Chuine I, Kirkpatrick M, Ronce O (2012) How do genetic correlations affect species range shifts in a changing climate? Ecol Lett 15(3):251–259. 22. Kubisch A, Degen T, Hovestadt T, Poethke HJ (2013) Predicting range shifts under global change: The balance between local adaptation and dispersal. Ecography 36(8): 873–882. 23. García-Ramos G, Kirkpatrick M (1997) Genetic models of rapid evolutionary divergence in peripheral populations. Evolution 51(1):21–28. 24. Kirkpatrick M, Barton NH (1997) Evolution of a species’ range. Am Nat 150(1):1–23. 25. Schiffers K, Bourne EC, Lavergne S, Thuiller W, Travis JMJ (2013) Limited evolutionary rescue of locally adapted populations facing climate change. Philos Trans R Soc Lond B Biol Sci 368(1610):20120083. 26. Lopez S, Rousset F, Shaw FH, Shaw RG, Ronce O (2008) Migration load in plants: Role of pollen and seed dispersal in heterogeneous landscapes. J Evol Biol 21(1):294–309. 27. Ouborg NJ, Piquot Y, Groenendael JMV (1999) Population genetics, molecular markers and the study of dispersal in plants. J Ecol 87(4):551–568. 28. Petit RJ, et al. (2005) Comparative organization of chloroplast, mitochondrial and nuclear diversity in plant populations. Mol Ecol 14(3):689–701.

E5748 | www.pnas.org/cgi/doi/10.1073/pnas.1607612113

29. Savolainen O, Pyhäjärvi T, Knürr T (2007) Gene flow and local adaptation in trees. Annu Rev Ecol Evol Syst 38(1):595–619. 30. Nagylaki T (1975) Conditions for the existence of clines. Genetics 80(3):595–615. 31. Barton NH (2001) Adaptation at the edge of a species’ range. Integrating Ecology and Evolution in a Spatial Context, eds Silvertown J, Antonovics J (Blackwell Science, New York), pp 365–392. 32. Mirrahimi S, Raoul G (2013) Dynamics of sexual populations structured by a space variable and a phenotypical trait. Theor Popul Biol 84(1):87–103. 33. Johnson T, Barton N (2005) Theoretical models of selection and mutation on quantitative traits. Philos Trans R Soc Lond B Biol Sci 360(1459):1411–1425. 34. Crawford TJ (1984) The estimation of neighbourhood parameters for plant populations. Heredity 52(2):273–283. 35. Maynard Smith J (1968) “Haldane’s dilemma” and the rate of evolution. Nature 219(5159):1114–1116. 36. Bulmer MG (1980) The Mathematical Theory of Quantitative Genetics (Clarendon Press, Oxford). 37. Peterson AT, et al. (2011) Ecological Niches and Geographic Distributions (Princeton Univ Press, Princeton). 38. Bertrand R, Perez V, Gégout JC (2012) Disregarding the edaphic dimension in species distribution models leads to the omission of crucial spatial information under climate change: The case of Quercus pubescens in France. Glob Change Biol 18(8):2648–2660. 39. Wisz MS, et al. (2013) The role of biotic interactions in shaping distributions and realised assemblages of species: Implications for distribution modeling. Biol Rev Camb Philos Soc 88(1):15–30. 40. Catullo RA, Ferrier S, Hoffmann AA (2015) Extending spatial modelling of climate change responses beyond the realized niche: Estimating, and accommodating, physiological limits and adaptive evolution. Glob Ecol Biogeogr 24(10):1192–1202. 41. Kawecki TJ (2003) Sex-biased dispersal and adaptation to marginal habitats. Am Nat 162(4):415–426. 42. Polechová J, Barton NH (2015) Limits to adaptation along environmental gradients. Proc Natl Acad Sci USA 112(20):6401–6406. 43. Tufto J (2000) Quantitative genetic models of the balance between migration and stabilizing selection. Genet Res 76(3):285–293. 44. Huisman J, Tufto J (2012) Comparison of non-gaussian quantitative genetic models for migration and stabilizing selection. Evolution 66(11):3444–3461. 45. Bocedi G, et al. (2013) Effects of local adaptation and interspecific competition on species’ responses to climate change. Ann N Y Acad Sci 1297(1):83–97. 46. Louthan AM, Doak DF, Angert AL (2015) Where and when do species interactions set range limits? Trends Ecol Evol 30(12):780–792. 47. Hargreaves AL, Bailey SF, Laird RA (2015) Fitness declines toward range limits and local adaptation to climate affect dispersal evolution during climate-induced range shifts. J Evol Biol 28(8):1489–1501. 48. Hoffmann AA, Sgrò CM (2011) Climate change and evolutionary adaptation. Nature 470(7335):479–485. 49. Galassi M, et al. (2009) GNU Scientific Library Reference Manual (Network Theory Ltd., Surrey, UK), 3rd Ed.

Aguilée et al.

Supporting Information Aguilée et al. 10.1073/pnas.1607612113 Appendix S1. Numerical Resolution of Eqs. 1 and 3 We used numerical resolution to evaluate the stability of the equilibria analytically found with the simplified model, and to determine the equilibria of the model with local density dependence and evaluate their stability. Using variables ζ = x − vt and ηðζ, tÞ = logðnðζ, tÞÞ, Eq. 1 becomes the following: ∂t η = v  ∂ζ η +

σ 2s 2




∂ζ,ζ η + ð∂ζ ηÞ2



  ~ ðz − bζÞ2 Vp n − − , + r0 1 − k 2Vs 2Vs

R ~ = expðηðζ, tÞÞdζ for the simplified model and n ~= where n expðηðζ, tÞÞ for the model with local density dependence. For both models, Eq. 3 becomes the following: ∂tz = v  ∂ζz +

z − bζ σ   ∂ζ,ζz + σ 2   ∂ζz  ∂ζ η − Vg . Vs 2 2

These equations can be numerically integrated using the finitedifferences method (forward differences at the extreme left of the domain, backward differences at the extreme right of the domain, and centered finite differences elsewhere). We used the explicit embedded Runge–Kutta–Fehlberg (4, 5) method with adaptive step as implemented in the GNU Scientific Library (49). The domain of integration in space was defined large enough so that, if the solution is a traveling wave, it does not go out of the domain before the end of the integration. The step in space was δx = σ s =4. Numerical integration was stopped when the population reached an equilibrium defined as no measurable variation of total population size during a minimum of 100 consecutive time steps of integration. When the equilibrium reached was a traveling wave, we continued the numerical integration over 50 time units to measure the speed of the spatial density shift and of the ecological niche shift. We tested a wide range of initial conditions. For the initial phenotypic cline, we used the following: a flat phenotypic cline [zðx, 0Þ = 0], a phenotypic cline with a slope close to the slope predicted by the simplified model, a phenotypic cline shifted to the left or to the right compared with the analytical prediction from the simplified model, a phenotypic cline corresponding the prediction of the simplified model with white noise added. For the initial population size, we used the following: the same population size on the whole domain of integration, one narrow pic of population density at ζ close to 0, a population size as predicted by the simplified model but shifted to the left or to the right of the domain, and a population size corresponding to the prediction of the simplified model with white noise added. When the equilibrium reached was a traveling wave, the slope s of the phenotypic cline was estimated using a linear regression pffiffiffiffiffi weighted by local population size. The size of the range vn was estimated as the SD of the geographical position of individuals. The rate of the spatial density shift cn was estimated as the speed of the position of the maximum of population density (which may not coincide with the mean position of individuals, because numerical resolution allows for nonsymmetric solutions). The spatial lag of the phenotypic cline −ln was estimated as the distance between the maximum of population density and the location where the mean phenotype at the maximum of population density is optimal. The rate of ecological niche shift d was estimated as the speed at which the mean phenotype at the maximum of population density changes. Aguilée et al. www.pnas.org/cgi/content/short/1607612113

Appendix S2. Model with Evolving Genetic Variance We consider here a population structured by a genotypic value g and a space variable x. We denote ηðx, t, gÞ the density of individuals with a genotypic value g at time t at location x. Following Mirrahimi and Raoul (32), the evolution of population size is described by the following: ∂t ηðx, t, gÞ =

σ 2s ∂x,x ηðx, t, gÞ 2

nðx, tÞ + − ηðx, t, gÞ ðβ − r0 Þ + r0 k

Z

ðz − θðx, tÞÞ2 fg ðzÞdz 2Vs

!

+ Sðx, t, gÞ, [S1] where nðx, tÞ is the local population density, β is the per capita birth rate, fg ðzÞ is the distribution of phenotypes among individuals with genotype g, and Sðx, t, gÞ is the number of seeds with genotype g produced at location x at time t (notation is otherwise similar to the one used in the other models; Table 1). The first term of Eq. S1 describes the demographic effect of dispersal of seeds only; we assume that pollen does not limit population growth. The second term of Eq. S1 describes death: (i) natural death and death due to density dependence, which is assumed to be local [the per capita growth rate under a logistic model is r0 ð1 − nðx, tÞ=kÞ = β − ðβ − r0 Þ − r0 nðx, tÞ=k; in Eq. S1, ðβ − r0 Þ is thus the natural death rate, and r0 nðx, tÞ=k is the death rate due to density dependence], and (ii) death due to the mismatch between the individual phenotype and the local optimal phenotype (i.e., natural selection; integral term in Eq. S1). We assume that the distribution fg ðzÞ of phenotypes among individuals with genotype g is Gaussian with pffiffiffiffiffiffiffiffiffiffi ffi mean g and variance Ve: fg ðzÞ = expð−ðz − gÞ2 =ð2Ve ÞÞ= 2πVe. As mentioned above, the third term of Eq. S1, Sðx, t, gÞ, is the number of seeds with genotype g produced at location x at time t: it describes reproduction. At time t, the reproduction term can be decomposed as Sðx, t, gÞ = βnðx, tÞPðx, gÞ where βnðx, tÞ is the amount of ovules at location x and Pðx, gÞ is the probability that an ovule at location x gives a seed with genotype g. Conditioning by the genotype of the parents, this term can be decomposed as follows: ZZ Pðx, gÞ = Qðg, g\ , g_ ÞP mother ðx, g\ ÞP father ðx, g_ Þdg\ dg_ , where Qðg, g\ , g_ Þ is the probability that an ovule at location x gives a seed with genotype g knowing that the mother has genotype g\ and the father has genotype g_, P mother ðx, g\ Þ is the probability that the mother at location x has genotype g\, and P father ðx, g_ Þ is the probability that the father of the seed to be born at location x has genotype g_. We use here the infinitesimal model (36) to describe the genetic transmission of the trait: phenotype depends on a large number of loci with infinitesimal and additive effects, which implies that allele frequency changes can be neglected, so that the genetic variance at linkage equilibrium, vLE, remains constant. The genotype of a newborn individual is distributed according to a Gaussian kernel centered on the mean genotype of the parents and with variance pffiffiffiffiffiffiffiffiffiffiffiffivLE, that is, Qðg, g\ , g_ Þ = expð−ðg − ðg\ + g_ Þ=2Þ2 =ð2vLE ÞÞ= 2πvLE. This model of inheritance approximates well more realistic explicit multilocus models (43, 44). The term P mother ðx, g\ Þ is the proportion of 1 of 8

individuals with genotype g\ at location x: P mother ðx, g\ Þ = ηðx, t, g\ Þ=nðx, tÞ. The term P father ðx, g_ Þ is the proportion of pollen from a father with genotype g_ available at location x. We denote Pðx, x′Þ the probability that an individual at location x′ sends pollen at location x. Assuming a Gaussian dispersal kernel with variance σ 2p for pollen, Pðx, x′Þ = expð−ðx − x′Þ2 =ð2σ 2p ÞÞ= qffiffiffiffiffiffiffiffiffiffi R 2πσ 2p. Then, we can write P father ðx, g_ Þ = ηðx′, t, g_ ÞPðx, x′Þdx′= R nðx′, tÞPðx, x′Þdx′. Finally, computing the integral in Eq. S1 and simplifying the expression of Sðx, t, gÞ, Eq. S1 becomes the following: σ 2s ∂x,x ηðx, t, gÞ 2   nðx, tÞ ðg − θðx, tÞÞ2 Ve − − + ηðx, t, gÞ −ðβ − r0 Þ − r0 k 2Vs 2Vs R ZZ ηðx′, t, g_ ÞPðx, x′Þdx′ dg\dg_. +β Qðg, g\ , g_ Þηðx, t, g\ Þ R nðx′, tÞPðx, x′Þdx′ [S2]

∂t ηðx, t, gÞ =

From Eq. S2, we can retrieve the R macroscopic variables: the local population size is nðx, tÞ = ηðx, t, gÞdg, the mean local phenotypic R trait zðx, tÞ equals the mean local genotypic value gðx, tÞ = gηðx, t, gÞ=nðx, tÞdg, and the local additive genetic variR ance is Vg ðx, tÞ = ðg − gðx, tÞÞ2 ηðx, t, gÞ=nðx, tÞdg. In this model, the variance at linkage equilibrium vLE is constant in space and in time; the additive genetic variance can nevertheless change due to changes in linkage disequilibrium, and the full distribution of offspring (i.e., unconditional to the parental genotypes) is not constrained to be Gaussian. We solved Eq. S2 numerically using the finite-differences method (forward differences at the extreme left of the domain, backward differences at the extreme right of the domain, and centered finite differences elsewhere) with the Euler method. The domain of integration in space and in genotype was defined large enough so that, if the solution is a traveling wave, it does not go out of the domain before the end of the integration. The pffiffiffiffiffiffiffi step in space was δx = σ s =4, the step in genotype was δg = vLE =2, and the step in time was δt = 0.1δx. Numerical resolution was stopped when the population reached an equilibrium defined as no measurable variation of total population size during a minimum of 100 consecutive time steps of integration. Numerical resolution was initialized with a Gaussian distribution of individuals in space with variance σ 2s. The distribution of genotypes was also Gaussian, with variance 2vLE, and centered on the optimal phenotype at the center of the geographical range. When the equilibrium reached was a traveling wave, the parameters s, vn, cn, ln, and d were estimated as described in Appendix S1. Numerical Resolution of Eqs. 1 and 3. Appendix S3. Resolution of the Simplified Model: TravelingWave Equilibrium S3.1. Expression of the Parameters of the Traveling-Wave Equilibrium.

Let us assume that there is a traveling-wave solution for Eq. 4, which can be written as follows: ! ðX − Cn T − Ln Þ2 NðX, TÞ = N0 exp − 2Vn , [S3] ZðX, TÞ = SðX − Cz T − Ln Þ. This solution is equivalent to the one given in the main text (Eq. 6) with Cz = Cn − D=S. The solution must satisfy Eq. 4 for all values Aguilée et al. www.pnas.org/cgi/content/short/1607612113

of X and T, which leads, after substituting Eq. S3 into Eq. 4, to the following: 2S − AðS − BÞ = 0, ð1 − γÞVn

[S4]

2SCn + ASCz − ABV = 0, ð1 − γÞVn

[S5]

Cz +

1−Λ−

2Ln + ALn = 0, ð1 − γÞVn

  ðSLn Þ2 1 Ln Ln − + Cn + = 0, Vn Vn 2 Vn

[S6]

[S7]

  Cn Ln Cn + 2 + SLn ðBS − SCz Þ = 0, Vn Vn

[S8]

 2 Cn ðBV − SCz Þ − 2 = 0, Vn

[S9]

  1 Ln Cn + 2 = 0, SLn ðS − BÞ − Vn Vn

[S10]

2 − ðS − BÞ2 = 0, Vn2

[S11]

2

ðS − BÞðBV − SCz Þ + 2

Cn = 0. Vn

[S12]

Parameters for which we look for an expression are S, Vn, Cz, Cn, Ln, and Λ. Let us assume for now that S − B > 0. Then, from Eq. S11: S−B=

pffiffiffi 2 . Vn

[S13]

Using this expression and Eq. S4 gives the following: A S = −pffiffiffi ð1 − γÞ. 2

[S14]

Combining Eqs. S13 and S14 leads to the following: Vn =

2 pffiffiffi . −B 2 − Að1 − γÞ

[S15]

Using Eq. S14, Eq. S6 can be rewritten as follows: ð1 − γÞCz + 2

pffiffiffi Ln + SLn 2 = 0. Vn

[S16]

Using Eq. S13, Eq. S10 can be rewritten as follows: pffiffiffi Ln SLn 2 + Cn + 2 = 0. Vn

[S17]

Subtracting Eqs. S16 and S17 leads to the following: Cz =

Cn . 1−γ

[S18]

2 of 8

Substituting Eqs. S14, S15, and S18 into Eq. S5 gives after simplification the following: Cn =

V . 1 + −BApffiffi2 γ

[S19]

Combining Eqs. S18 and S19 leads to the following: Cz =




V

ð1 − γÞ 1 + −BApffiffi2 γ

.

[S20]

From Eq. S6, into which we substitute Eqs. S15 and S20, we get the following: V . Ln = − pffiffiffi −B 2 + Aγ

[S21]

Using Eqs. S18 and S15 into Eq. S6 gives the following: pffiffiffi Cn = BLn 2. [S22] Finally, pwe ffiffi can obtain Λ from Eq. S7: using Eq. S22, then that S = B + Vn2 (Eq. S13), Λ can be written as follows: Λ=1−

1 ðBLn Þ2 , − 2 Vn

[S23]

where Vn and Ln are defined by Eqs. S15 and S21. We checked that Eqs. S4–S12 are all verified using the above expression of S, Vn, Cn, Cz, Ln, and Λ. Repeating the resolution assuming that S − B < 0 leads to the expression of S, Vn, Cn, and Ln given in the main text and leads to the following expression for Cz: Cz =




V

ð1 − γÞ 1 + jBjApffiffi2 γ

.

[S24]

The expression of D in the main text can then be computed using D = SðCn − Cz Þ (Eq. 6 and Eq. S3). S3.2. Mean Growth Rate. The center of the range, where the size of the population is maximal, is located in X1 = Cn T + Ln (Eq. 6). Using Eq. 5 and Eq. S3, we find that the mean growth rate R reaches a maximum at location X2 = ðTðBV − Cz SÞ − SLn Þ=ðB − SÞ. Then, substituting S, Cn, Ln, and Cz with their expressions (Eqs. 7, 9, and 11, and Eq. S24), the distance between the maximum of the mean growth rate and the p maximum ffiffiffi pof ffiffiffi population density pffiffiffi can be written as X2 − X1 = V jBj 2=ððjBj 2 − ð1 − γÞAÞðjBj 2 + AγÞÞ. Assuming that the environmental gradient shifts in space (V > 0), andpbecause at the traveling-wave equilibrium we necessarily have ffiffiffi jBj 2 − ð1 − γÞA > 0 (Eq. 14; Appendix S3. Resolution of the Simplified Model: Traveling-Wave Equilibrium, S3.4. Critical Thresholds), we have X2 > X1, that is, the maximum of population size lags behind the maximum of the mean growth rate. Furthermore, the above expression of X2 − X1 shows that it is a decreasing function of γ, that is, the highest mean growth rate is closer to the center of the range when pffiffiffi pollen pffiffiffiis dispersed further. Note also that X2 − X1 = −Ln jBj 2= ðjBj 2 − ð1 − γÞAÞ with Ln given by Eq. 11, that is, the spatial lag of the maximum of population density behind the maximum of the mean growth rate is proportional to the spatial lag of the phenotypic cline Ln. S3.3. Spatial Lag of the Phenotypic Cline. Individuals at the center of

the range are located in X1 = Cn T + Ln and have a mean phenotype ZðX1 , TÞ = SðCn − Cz ÞT (Eq. S3). Let X2 be the location Aguilée et al. www.pnas.org/cgi/content/short/1607612113

where this mean phenotype ZðX1 , TÞ would be perfectly adapted. Then X2 is such that SðCn − Cz ÞT = BðX2 − VTÞ. Using Eqs. 7 and 9 and Eq. S24, we can check that SðCn − Cz Þ = BðCn − V Þ. Thus, BðCn − V ÞT = BðX2 − VTÞ, which allows to write X2 = Cn T. Finally, we get X2 − X1 = −Ln. Therefore, −Ln is the constant spatial distance that the population would have to move for individuals at the center of the range to be perfectly adapted. S3.4. Critical Thresholds. The traveling-wave solution may exist only if (i) the total population size is positive, that is, Λ > 0, and if (ii) the size of the range is not null, that is, Vn > 0. Using Eq. S23 and the expressions of Vn and Ln (Eqs. 8 and 11), the first condition can crit = be written as a critical rate of environmental change: V < VTW ffi pffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffi 2ð1 + Aγ=ðjBj 2 ÞÞ 1 − jBj 2=2 + Að1 − γÞ=2 (Eq. 12). The two occurrences of γ in this critical rate have antagonistic effects. By crit manually computing VTW , we can track their origin. The first occrit , comes from the currence of γ, which has a positive effect on VTW term Ln in Eq. S23. This leads us to conclude that increasing the distance of pollen dispersal may allow to track a faster environmental shift because dispersing pollen further reduces the spatial lag of the phenotypic cline (Eq. 11). The second occurrence of γ, crit , comes from the term Vn in Eq. which has a negative effect on VTW S23. We thus conclude that the ability to track a faster environmental shift may be compromised by dispersing pollen further because this contracts the range of the population (Eq. 8). The condition (ii) for the existence of the traveling-wave equilibrium can be written using Eq. 8 either as a critical pffiffiffi slope of 14] the environmental gradient [jBj > jBcrit ffiffiffi or TW j = Að1 − γÞ= 2; Eq.p as a critical distance of pollen dispersal p (γffiffiffi> γ crit TW = 1 − jBj 2=A; Eq. 15). Given that S = signðBÞAð1 − γÞ= 2 (Eq. 7), jSj = jBcrit TW j. Therefore, the phenotypic cline is necessarily flatter than the environmental gradient at the traveling-wave equilibrium (jSj < jBj). This condition [also found in previous studies without pollen dispersal (19, 24)] means that the genetic swamping effect of dispersal must be sufficient to prevent the range to expand without limit (i.e., the uniform-density equilibrium). opt crit of γ (Eq. 12) has a maximum at γ = γ TW = The function VTW pffiffiffi 2=3 − ðjBj 2 − 4=3Þ=A (Eq. 13). Biologically possible values of γ crit are those between 0 and 1; we consider the variations of VTW only on this interval. Consequently, the critical rate of environcrit always increases with the distance of pollen mental change VTW opt dispersal γ if γ TW ≥ 1, which can be rewritten using Eq. 13 as pffiffiffi jBj ≤ ð4 − AÞ=ð3 2 Þ. In other words, when the environmental gradient is shallow, the optimal distance of pollen dispersal is as large as biologically possible (γ → 1). Conversely, the optimal contribution of opt pollen dispersal is no dispersal (γ = 0) when γ TWp≤ ffiffiffi 0, which corresponds to a steep environmental gradient [jBj ≥ 2ð2 + AÞ=3]. The optimal distance of pollen dispersal is intermediate when the slope of the environmental gradient is intermediate.

Appendix S4. Traveling-Wave Equilibrium: Other Parameters Effect Here, we detail the effect of the genetic potential for adaptation A, of the slope B of the environmental gradient, of the rate V of environmental change, and of the distance of seed dispersal σ s on the traveling-wave equilibrium. The effects described below can be deduced from Eqs. 7–15, and from the corresponding equations using nonscaled variables (Table S1) for the effect of the distance of seed dispersal σ s. The results described below are consistent with those found in previous analytical studies without pollen dispersal (Appendix S5. Traveling-Wave Equilibrium: Comparison with Duputié et al. (21)). Slope of the Phenotypic Cline S. Independently of the distance of pollen dispersal, the phenotypic cline is steeper as the genetic potential for adaptation A increases, but it depends neither on the value of the slope B of the environmental gradient (except for its sign) nor 3 of 8

on the shifting rate of the environmental gradient V. A higher distance of seed dispersal σ s flattens the realized phenotypic cline. Size of the Spatial Range Measured by Vn. The size of the spatial range does not depend on the speed of environmental change V. A steeper environmental gradient (higher jBj) contracts the spatial range, whereas a better potential for adaptation (higher A) expands it. The size of the spatial range reaches a minimum for an intermediate distance of seed dispersal σ s. Shifting Rate of the Population Density Cn. The spatial shift of the population density occurs at a rate proportional to the rate of environmental shift V. When pollen is dispersed (γ > 0), it is slowed by a better adaptive potential A and accelerated by a steeper environmental gradient (higher jBj) and by a larger distance of seed dispersal σ s. Rate of Ecological Niche Shift D. For any non-null distance of pollen

dispersal (γ > 0), the ecological niche shift occurs at a rate proportional to the rate of environmental shift V. It is faster when either the genetic potential for adaptation A increases or when the slope jBj of the environmental gradient increases. The rate of ecological niche shift is slower when seeds disperse further (higher σ s). Spatial Lag of the Phenotypic Cline Ln. The spatial lag of the phenotypic cline is proportional to the rate of environmental shift V. The phenotypic cline tracks the environmental gradient more closely when the environmental gradient is steeper (higher jBj), when seeds are dispersed further (higher σ s), and, assuming that the distance of pollen dispersal is not null (γ > 0), when the adaptive potential A is greater. Critical Rate of Environmental Change V crit TW . Independently of the

distance of pollen dispersal, the critical rate of environmental change is higher when the genetic potential for adaptation A increases and when the environmental gradient is shallower (lower jBj). The effect of the distance of seed dispersal σ s on the critical rate of environmental is nonmonotonous: there is an intermediate distance of seed dispersal that maximizes the sustainable rate of environmental change. Optimal Distance of Pollen Dispersal γopt TW . The optimal distance of pollen dispersal relative to the distance of seed dispersal increases for a higher genetic potential for adaptation A, for a lower slope jBj of the environmental gradient, and for a shorter distance of seed dispersal σ s. Critical Slope of the Environmental Gradient Bcrit TW and Critical Distance of Pollen Dispersal γcrit TW . A higher genetic potential for adaptation

A and a shorter distance of seed dispersal σ s increase both thresholds. The critical distance of pollen dispersal becomes higher when the environmental gradient is shallower (lower jBj).

Appendix S5. Traveling-Wave Equilibrium: Comparison with Duputié et al. (21) Here, we show that all our results of the simplified model converge to those in the study by Duputié et al. (21) when pollen does not disperse. In Duputié et al. (21), there is no pollen dispersal (γ = 0), the slope of the environmental gradient is positive (b > 0), and the weak selection assumption is used (Vg =Vs is small). Using nonscaled variables (Table S1), the slope of the phenotypic cline (Eq. 7) becomes the following: Vg s = signðbÞ pffiffiffiffiffi ð1 − γÞ. σ s Vs pffiffiffiffiffi Without pollen dispersal (γ = 0), s = Vg =ðσ s Vs Þ, which is the same expression as in the study by Duputié et al. (21) when we assume Aguilée et al. www.pnas.org/cgi/content/short/1607612113

b > 0. At any fixed distance of pollen dispersal γ, all other parameters have the same effect on s here and in the study by Duputié et al. (21). Using nonscaled variables, Eq. 8 becomes the following: vn =

σ2V pffiffiffiffiffi s s . σ s jbj Vs − Vg ð1 − γÞ

[S25]

2 Without pffiffiffiffiffi pollen dispersal and assuming b > 0, vn = σ s = ðσ s b= Vs − Vg =Vs Þ. Using Taylor series approximations to the first degree of Vg =Vs around 0 (i.e., using the weak selection assumption), we get the following: pffiffiffiffiffi   Vg Vs σ s + pffiffiffiffiffi , vn ≈ b b Vs

which is the same expression as in the study by Duputié et al. (21). With or without the weak selection assumption, vn increases when Vg increases and it decreases when b increases. However, with the weak selection assumption, the size of the range always increases with the distance of seed dispersal σ s, whereas with the exact expression (Eq. S25) vn is a nonmonotonous function of σ s: pffiffiffiffiffi it reaches a minimum at σ s = σes = 2Vg =ðb Vs Þ. In other words, when the genetic potential for adaptation (Vg) is high and the distance of seed dispersal is short, the size of the range decreases as the distance of seed dispersal increases. Using nonscaled variables, Eq. 9 becomes the following: cn =

v g ffiffiffiffi γ 1 + σ jbjp V

V

s

,

s

and Eq. 10 becomes the following: bvVg γ pffiffiffiffiffi . d=− σ s jbj Vs + Vg γ Without pollen dispersal, cn = v and d = 0, which are the same expressions as in the study by Duputié et al. (21). When pollen disperses (γ > 0), both rates are proportional to the rate of environmental shift v; they are both accelerated by a steeper environmental gradient (higher jbj) and/or a higher distance of seed dispersal σ s; a better adaptive potential (Vg) slows the spatial range shift cn and accelerates the ecological niche shift d. The lag of the phenotypic cline (Ln; Eq. 11) is not defined similarly here and in the study by Duputié et al. (21). Here, −Ln is the spatial distance that the population would have to move for individuals at the center of the range to be perfectly adapted (Appendix S3. Resolution of the Simplified Model: Traveling-Wave Equilibrium, S3.3. Spatial Lag of the Phenotypic Cline); in the study by Duputié et al. (21), −Ln is the spatial distance between the center of the range and the location where individuals are currently perfectly adapted. Let us denote −L′n this spatial distance between the center of the range and the location where individuals are currently perfectly adapted in our model. It is computed in Appendix S3. Resolution of the Simplified Model: Traveling-Wave Equilibrium, S3.2. Mean Growth Rate as follows: pffiffiffi V jBj 2 L′n = − pffiffiffi  pffiffiffi . Bj 2 − ð1 − γÞA jBj 2 + Aγ This expression can be written using nonscaled variables as follows: v l′n = − bσ Vg , pffiffiffisffi − V Vs

s

when pollen does not disperse and assuming b > 0. Using Taylor series approximations to the first degree of Vg =Vs around 0, we get the following: 4 of 8

  v pffiffiffiffiffi Vg Vs + , l′n ≈ − bσ s bσ s which is the same expression as ln in the study by Duputié et al. (21). Using nonscaled variables, the critical rate of the environmental shift above which the traveling-wave equilibrium does not exist (Eq. 12) becomes the following: vcrit TW

 = σs +

Vg pffiffiffiffiffi γ jbj Vs

ffi  sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  Vp jbjσ s Vg − pffiffiffiffiffi + ð1 − γÞ. 2 r0 − 2Vs Vs Vs

Without pollen dispersal and assuming b > 0, vcrit TW = qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffi σ s 2ðr0 − Vp =ð2Vs ÞÞ − bσ s = Vs + Vg =Vs, which is the same expression as in the study by Duputié et al. (21). At any fixed distance of pollen dispersal γ, all other parameters have the same effect on vcrit TW here and in the study by Duputié et al. (21). Appendix S6. Resolution of the Simplified Model: Uniform-Density Equilibrium Let us assume that there is a solution for Eq. 4, which can be written as follows: NðX, TÞ = N0 , ZðX, TÞ = SðX − VTÞ + Lz .

[S26]

This solution must satisfy Eq. 4 for all values of ðX − VTÞ, which leads, after substituting Eq. S26 into Eq. 4, to the following: ðS − BÞ2 = 0,

[S27]

Lz ðS − BÞ = 0,

[S28]

L2z − 2ð1 − ΛÞ = 0,

[S29]

AðS − BÞ = 0,

[S30]

ALz − SV = 0.

[S31]

From Eq. S30, S = B (Eqs. S27 and S28 are then verified); from Eq. 2 S31, Lz = BV A ; and from Eq. S29, Λ = 1 − ð1=2ÞðBV =AÞ . The uniform-density equilibrium exists when the scaled total population pffiffiffi crit size Λ is positive, that is, from Eq. S29, when V < VUD = A 2=jBj. At the uniform-density equilibrium, the mean phenotypic cline has the same slope B as the environmental gradient and moves across space at the same rate V as the environmental gradient. Maladaptation is constant in space and time and equals jBjV =A. None of the formulas above depends on γ: the properties of the uniform-density equilibrium are consequently independent of the distance of pollen dispersal. Appendix S7. Robustness of the Results S7.1. Results with Local Density Dependence. To understand the effect of local density dependence on the properties of the traveling-wave equilibrium, let us first accurately define the center of the spatial range of the population: the range is the part of the cumulative distribution function of the population density between quantiles 0.025 and 0.975 (i.e., the central part where 95% of the individuals are found); its center is here defined by the average of these quantiles. Under local density dependence, competition reduces less the local mean growth rate at the margins, where local density is the lowest, than at the center of the range, where local density is higher (Eq. 2). Consequently, compared with global density deAguilée et al. www.pnas.org/cgi/content/short/1607612113

pendence for which the maximum of population density is at the center of the range and considering locations where the population will not decline in the future (i.e., the leading edge, where the mean growth rate is positive; Fig. 1), local density dependence shifts the maximum of population density in the direction of the environmental change (Fig. S1 A and C). The resulting asymmetrical spatial distribution of individuals induces a swamping effect of dispersal stronger at the leading edge and a swamping effect weaker at the trailing edge because the gradient in density is more shallow. The slope of the realized phenotypic cline, measured as the coefficient of linear regression weighted by local population density, is therefore closer to the slope of the environmental gradient when density dependence is local than when it is global (Fig. S1 B and D and Fig. 5A). Note that, as there is a spatial lag of the phenotypic cline (Ln ≠ 0, Eq. 11), the stronger swamping effect at the leading edge that flattens the realized phenotypic cline at the leading edge compared with the center of the range makes the phenotype of individuals at the leading edge more similar to the optimal phenotype (Fig. S1D). The shallower realized phenotypic cline under local density dependence has cascading effects on the features of the traveling wave (Fig. 5 B–D): the traveling wave is wider than expected from the simplified model, the population density shifts in space at a higher rate [the population density may shift in space faster than the environmental gradient when the distance of pollen dispersal is not high, as found by Polechová et al. (19) without pollen dispersal], and the lag of the phenotypic cline is reduced. With local density dependence, the ecological niche shift is pulled in a different direction; the resulting shift will be slower or faster in magnitude depending on what the simplified model predicted (Fig. 5E). The quantitative effect of local density dependence is weak for crit parameter values close to the extinction thresholds VTW and γ crit TW (Fig. S1 E and F). Indeed, when the population is close to extinction, local density dependence is close to zero over the whole range, and competition becomes homogeneous in space. The model with local density dependence then converges to the simplified model for which the spatial distribution of individuals at equilibrium is a Gaussian traveling wave. S7.2. Results with Evolving Genetic Variance. Allowing the additive genetic variance to change due to changes in linkage disequilibrium (Materials and Methods, Model with Evolving Genetic Variance), we observed that it indeed changes through space and time: when the population reaches the traveling-wave equilibrium, the spatial profile of the genetic variance reaches an equilibrium with a strong decrease at both margins of the population (Fig. S2 A and C); when the population tends to the uniform-density equilibrium, the genetic variance increases over time on the whole space during the course of numerical resolution (Fig. S2 B and D). Despite these spatial and temporal variations of the genetic variance, the existence and the local stability of both equilibria are rather similar to the results found with the model with local density dependence and fixed genetic variance (Fig. 4C). Nevertheless, we found no parameter values for which both the traveling-wave equilibrium and the uniform-density equilibrium existed simultaneously (Fig. 4C). When the population is at the traveling-wave equilibrium, we measured the additive genetic variance Vg as the mean genetic variance over space weighted by local population size. The distance of pollen dispersal may have opposite effects on the genetic variance because, on the one hand, it brings alleles from populations differentiated along the realized phenotypic cline, which may increase the local genetic variance, but on the other hand, pollen dispersal flattens the phenotypic cline (Eq. 7) so that the genetic variance between populations is reduced. For all of the parameter values we investigated (Table 1), we found that the genetic variance Vg slightly decreased when the distance of pollen dispersal increased (Fig. S3). The second effect of pollen dispersal seems thus to dominate, at least under the infinitesimal model of trait inheritance and for the range of numerical values corresponding to bud set date in Picea sitchensis. 5 of 8

Fig. S1. Population density n (A, C, and E) and mean phenotype z (solid lines in B, D, and F) as a function of space ζ = x − vt with the simplified model (A and B, v = 32 km·gen.−1) and with the model with local density dependence (C and D: v = 32 km·gen.−1; E and F: v = 45 km·gen.−1) at the traveling-wave equilibrium. The dashed line shows the environmental gradient. The dotted lines show the linear regression weighted by local population size used to estimate the slope s of the realized phenotypic cline. The leading edge corresponds to the right side of each panel. Parameter values: default parameter values (Table 1), and r0 = 2 gen.−1, Vp = 2.5 d2, Vg = 2 d2, and γ = 0.85. We use here nonscaled parameters (Table S1). The third row (E and F) corresponds to a case close to extinction crit = 46.2 km·gen.−1), and the other rows (A and B, and C and D) show a case far from extinction. (vTW

Aguilée et al. www.pnas.org/cgi/content/short/1607612113

6 of 8

Fig. S2. Spatial profile at three different times (A and B) and dynamics (C and D) of the genetic variance with the model with evolving genetic variance when the population reaches the traveling-wave equilibrium (A and C) or tends to the uniform-density equilibrium (B and D). The mean genetic variance in C and D is the local genetic variance averaged over space weighted by local population size. Parameter values: default parameter values (Table 1), and r0 = 2 gen.−1, vLE = 1 d2. Traveling-wave equilibrium: v = 28.4 km·gen.−1, γ = 0.95; uniform-density equilibrium: v = 14.2 km·gen.−1, γ = 0.6.

Fig. S3. Genetic variance at the traveling-wave equilibrium (model with evolving genetic variance) as a function of the distance of pollen dispersal γ for three rates of the environmental shift. Parameter values: default parameter values (Table 1), and r0 = 1.5 gen.−1, vLE = 0.05 d2.

Aguilée et al. www.pnas.org/cgi/content/short/1607612113

7 of 8

Table S1. Correspondence between nonscaled and scaled variables, parameters, and outputs Nonscaled

Scaled

r0

V R0 = r0 − 2Vps

x

X=

t

T = R0 t

k

K = k Rr00

n

N = Kn

 z

z ffi Z = pffiffiffiffiffiffiffi R V

λ

Λ = Kλ

Vg

A = R0 Vg s

b

ffi B = b R pσsffiffiffiffiffi 0 2Vs pffiffi V = v σ p2ffiffiffiffi R

pffiffiffiffiffiffi 2R0 σs

0

x

s

V

v 

r = r0 1 − nk~

s



− vtÞÞ2 − ðz − bðx 2Vs

V − 2Vps

0

R = 1 − Λ − 12ðZ − BðX − VT ÞÞ2


 2 n t − ln Þ nTW = n0 exp −ðx − c2v n


 2 n T − Ln Þ NTW = N0 exp −ðX − C2V n

 zTW = sðx − cn t − ln Þ + dt

ZTW = SðX − Cn T − Ln Þ + DT

ffi ð1 − γÞ s = signðbÞ σ pgffiffiffi V

S = signðBÞ pAffiffi2 ð1 − γÞ

V

s

s

σ2 V

s s vn = σ jbjpffiffiffiffi V −V s

cn =

s

v 1+

Vn = jBjpffiffi2 −2Að1 − γÞ

g ð1 − γÞ

Cn = 1 + VApffi γ

pffiffiffi γ

Vg

σ s jbj

jBj 2

Vs

bvVg γ ffi d = −σ jbjpffiffiffi Vs + Vg γ s

pffiffi D = −jBjABVγ 2 + Aγ

ffiffiffiffis ln = −σ jbjpvV V +V

Ln = −jBjpVffiffi2 + Aγ

s

s

gγ

ffi
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
 V V V crit ffiγ pffiffiffiffis + g ð1 − γÞ vTW = σ s + jbjpgffiffiffi 2 r0 − 2Vps − jbjσ Vs V V s

s

ffi
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffi crit VTW = 2 1 + jBjApffiffi2 γ 1 − jBj2 2 + A2 ð1 − γÞ



 pffiffiffiffiffi Vp 1 4 γ opt + 23 TW = −Vg σ s jbj Vs − 3Vs r0 − 2Vs


pffiffiffi  1 4 2 γ opt TW = −A jBj 2 − 3 + 3

pgffiffiffiffi jbcrit TW j = σ V ð1 − γÞ

pAffiffi jBcrit TW j = 2 ð1 − γÞ

pffiffiffiffi σ s jbj Vs γ crit TW = 1 − Vg

pffiffi jBj 2 γ crit TW = 1 − A

nUD = n0

NUD = N0

s  zUD = bðx − vtÞ + bvV Vg

ZUD = BðX − VT Þ + BV A

V

s

Vg crit vUD =

s

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  Vp 2 r0 − 2Vs

pffiffiffiffi jbj Vs

pffiffi crit VUD = AjBj2

Population size (n, N) and mean phenotypic trait (z, Z) with a subscript TW (respectively UD) correspond to the travelingwave (respectively uniform-density) equilibrium.

Aguilée et al. www.pnas.org/cgi/content/short/1607612113

8 of 8