Is the sky the limit? On the expansion threshold of a species’ range

More than 100 years after Grigg’s influential analysis of species’ borders, the causes of limits to species’ ranges still represent a puzzle that has never been understood with clarity. The topic has become especially important recently as many scientists have become interested in the potential for species’ ranges to shift in response to climate change—and yet nearly all of those studies fail to recognise or incorporate evolutionary genetics in a way that relates to theoretical developments. I show that range margins can be understood based on just two measurable parameters: (i) the fitness cost of dispersal—a measure of environmental heterogeneity—and (ii) the strength of genetic drift, which reduces genetic diversity. Together, these two parameters define an ‘expansion threshold’: adaptation fails when genetic drift reduces genetic diversity below that required for adaptation to a heterogeneous environment. When the key parameters drop below this expansion threshold locally, a sharp range margin forms. When they drop below this threshold throughout the species’ range, adaptation collapses everywhere, resulting in either extinction or formation of a fragmented metapopulation. Because the effects of dispersal differ fundamentally with dimension, the second parameter—the strength of genetic drift—is qualitatively different compared to a linear habitat. In two-dimensional habitats, genetic drift becomes effectively independent of selection. It decreases with ‘neighbourhood size’—the number of individuals accessible by dispersal within one generation. Moreover, in contrast to earlier predictions, which neglected evolution of genetic variance and/or stochasticity in two dimensions, dispersal into small marginal populations aids adaptation. This is because the reduction of both genetic and demographic stochasticity has a stronger effect than the cost of dispersal through increased maladaptation. The expansion threshold thus provides a novel, theoretically justified, and testable prediction for formation of the range margin and collapse of the species’ range.


Introduction
Species' borders are not just determined by the limits of their ecological niche [1,2]. A species' edge is typically sharper than would be implied by continuous change in the species' environment (reviewed in [3, Table 2]). Moreover, although species' ranges are inherently dynamic, it is puzzling that they typically expand rather slowly [4]. The usual-but tautological-explanation is that lack of genetic variation at the range margin prevents further expansion [5]. Indeed, a species' range edge is often associated with lower neutral genetic variation [3,[6][7][8][9][10][11], suggesting that adaptive genetic variation may be depleted as well [12]. Yet why would selection for new variants near the edge of the range not increase adaptive genetic variance, thereby enabling it to continuously expand [5,13]? Haldane [14] proposed a general explanation: even if environmental conditions vary smoothly, 'swamping' by gene flow from central to marginal habitats will cause more severe maladaptation in marginal habitats, further reducing their population density. This would lead to a sharp edge to a species' range, even if genetic variance at the range margin is large. However, the consequences of dispersal and gene flow for evolution of a species' range continue to be debated [15][16][17][18]: a number of studies suggest that intermediate dispersal may be optimal [19][20][21][22][23]. Gene flow across heterogeneous environments can be beneficial because the increase of genetic variance allows the population to adapt in response to selection [13]. Current theory identifies that local population dynamics, dispersal, and evolution of niche-limiting traits (including their variance) and both genetic and demographic stochasticity are all important for species' range dynamics [13,[19][20][21][24][25][26][27][28]. Yet these core aspects have not been incorporated into a single study that would provide testable predictions for range limits in two-dimensional habitats.
As Haldane [14] previously pointed out, it is important to consider population and evolutionary dynamics across a species' range jointly, as their effects interact. Due to maladaptation, both the carrying capacity of the habitat and the population growth rate are likely to decrease-such selection is called 'hard' [29]. Classic deterministic theory [24] shows that when genetic variance is fixed, there are two stable regimes of adaptation to a spatially varying optimum (see Fig 1): (i) a 'limited adaptation', in which a population is only adapted to a single optimum or becomes a patchy conglomerate of discrete phenotypes, or (ii) continuous or 'uniform' adaptation, which is stable when the genetic variance, measured in terms of its cost in fitness (standing genetic load), is large relative to the maladaptation incurred by dispersal between environments (dispersal load). Under uniform adaptation, a species' range gradually expands-a stable boundary only forms when the genetic variance is too small to allow continuous adaptation to the spatially variable environment, and hence, limited adaptation is stable.
When genetic variance can evolve, such a limit no longer exists in infinitely large populations: the population maintains continuous adaptation as the environmental gradient steepens [13]. Deterministic theory thus predicts that a sharp and stable boundary to a species' range does not form when the environment changes smoothly. Uniform adaptation is the only stable regime when genetic variance can freely evolve in the absence of genetic drift [13], yet there is a limit to the steepness of the gradient. This limit arises because both the standing genetic load and the dispersal load increase as the gradient steepens, reducing the mean fitness (growth rate) of the population: when the mean fitness approaches zero, the population becomes extinct. Obviously, ignoring genetic drift is then unrealistic. In finite populations, genetic drift reduces local genetic variance [32], potentially qualitatively changing the dynamics. Indeed, it has been shown that for one-dimensional habitats (such as rivers), a sharp range margin arises when the fitness cost of dispersal across environments becomes too large relative to the efficacy of selection versus genetic drift [26]. However, most species live in two-dimensional habitats. There, allele frequencies can fluctuate over a local scale, as the correlations between them decline much faster across space than they do in linear habitats [33], and the effect of genetic drift changes qualitatively, becoming only weakly dependent on selection [34]. Is there still an intrinsic threshold to range expansion in finite populations when dispersal and gene flow occur over two-dimensional space rather than along a line? If so, what is its biological interpretation?

Results
I study the problem of intrinsic limits to adaptation in a two-dimensional habitat. Throughout, I assume that the species' niche is limited by stabilising selection on a composite phenotypic trait. This optimum varies across one dimension of the two-dimensional habitat-such as temperature and humidity with altitude. Demography and evolution are considered together. Selection is 'hard': both the rate of density-dependent population growth and the attainable equilibrium density decrease with increasing maladaptation. Both trait mean and genetic variance can freely evolve via change in allele frequencies and the associations among them (linkage disequilibria). The populations are finite, and both genetic and demographic stochasticity are included. The model is first outlined at a population level in terms of coupled stochastic differential equations. While it is not possible to obtain analytical solutions to this model, this formalisation allows us to identify the effective dimensionless parameters that describe the dynamics. Next, individual-based simulations are used to determine the driving relationship between the key parameters and test its robustness. The details are described in the Model section of the Methods.
The dynamics of the evolution of a species' range, as formalised by this model, are well described by three dimensionless parameters, which give a full description of the system. The first dimensionless parameter carries over from the phenotypic model [24]: the effective environmental gradient B measures the steepness of the environmental gradient in terms of maladaptation incurred by dispersal across a heterogeneous environment. The second parameter is the neighbourhood size of the population, N, which can be understood as the number of diploid individuals within one generation's dispersal range. Originally, neighbourhood size was defined by Wright [35] as the size of the single panmictic diploid population that would give the same probability of identity by descent in the previous generation. The inverse of neighbourhood size 1/ N hence describes the local increase of homozygosity due to genetic drift. The third dimensionless parameter is the ratio s/r Ã of the strength of selection s per locus relative to the strength of density dependence, r Ã . Detailed description of the parameters and their rescaling can be found in the Methods sections Parameters and Continuous model: Rescaling.
In order to see how the rescaled parameters capture the evolution of a species' range, I simulated 780 evolving populations, each based on different parameterisations, adapting to a linear gradient in the optimum. Depending on the parameters, the population either expands, gradually extending its phenotypic range by consecutive sweeps of loci advantageous at the edges, or the species' range contracts or disintegrates as adaptation fails. These two dimensionless parameters B and N give a clear separation between expanding populations, in which the neighbourhood size N is large relative to the effective environmental gradient B (shown in blue, Fig 2), and the rest, in which adaptation is failing. The separation gives an 'expansion threshold', estimated at N % 6.3B + 0.56 (red line). Above the expansion threshold, populations are predicted to expand (see Fig 3); below it, adaptation fails abruptly. If conditions change uniformly across space (as in these simulation runs, with constant gradient and carrying capacity), this means that adaptation fails everywhere-a species' range then either collapses from the margins (Fig 2, red hues) and/or disintegrates (Fig 2, open circles), forming a fragmented metapopulation (i.e., a spatially structured population consisting of discrete locally adapted subpopulations with limited dispersal among them).
When a metapopulation forms, it exhibits an extinction and colonisation dynamics. The subpopulations drift freely along the neutral spatial axis. Because the trait distributions of the subpopulations are unstable, the subpopulations also drift slowly along the environmental gradient. Over time, the metapopulation very slowly collapses to a virtually single trait value, with many subpopulaitons along the neutral axis. The subpopulations forming this metapopulation have only a very narrow phenotypic range and maintain locally only minimal adaptive variance. They correspond to the limited adaptation regime identified for a phenotypic model with genetic variance as a parameter [24]. In contrast to one-dimensional habitats [26], these  Fig 3). Open circles denote populations in which continuous adaptation has collapsed and the population consists of many discrete phenotypes adapted to a single optimum each (limited adaptation, Fig 4), whilst local genetic variance is very small. (Specifically, these are defined by mean heterozygosity smaller than 10% of the predicted value in the absence of genetic drift.) Simulations were run for 5,000 generations, starting from a population adapted to a linearly changing optimum in the central part of the available habitat. Populations that went extinct are marked with a black dot. Note that both axes are on a log scale. The top corner legend gives the colour-coding for the rate of range collapse and expansion in units of demes per generation; rates of collapse are capped at −1. The expansion threshold is fitted as a step function changing linearly along B: all blue dots are assigned a value of 1; all red dots and open circles are assigned a value of 0. The expansion threshold has a coefficient of determination R 2 = 0.94, calculated from 589 simulations (all but well-adapted stagnant populations). Data for this figure-and all subsequent ones-are deposited at Dryad Digital Repository, https://doi.org/10.5061/dryad.5vv37 [36].  closely matches the environmental gradient (grey) along the x-axis. (b) Population steadily expands, whilst population density stays continuous across space, with N " ¼ 19 AE 5:8 (mean ± standard deviation). The prediction at N ¼ 20 is shown by the blue contours; darker shading represents higher density. (c) Adaptation to the environmental gradient is maintained by a series of staggered clines: as each allele frequency changes from 0 to 1, the trait value increases by α. Population starts from the centre (blue hues reflect initial cline position relative to the centre of the range), and as it expands, new clines arising from loci previously fixed to 0 or 1 contribute to the adaptation (in red hues). At each location, multiple clines contribute to the trait (and variance); clines are shown at Y = 25. (d) Genetic variance changes continuously across space with mean V " G ¼ 0:032 AE 0:017 and stays slightly lower than is the deterministic prediction (green contours, V G = 0.045; higher variance is illustrated by darker shading). Deterministic predictions are based on [13] and are explained in the Methods section, along with the specification of the unscaled parameters. The population evolves for 2,000 generations, starting from a population adapted to the central habitat. The predicted neighbourhood size isN ¼ 34:6; effective environmental gradient is B = 0.48. https://doi.org/10.1371/journal.pbio.2005372.g003 Is the sky the limit? On the expansion threshold of a species' range Interestingly, the third dimensionless parameter s/r Ã has no detectable effect on the form of the expansion threshold. In other words, whilst the expansion threshold reflects the total fitness cost of dispersal in a heterogeneous environment, it appears independent of the strength of selection per locus s: the dashed lines in Fig 2 compare the estimated expansion threshold for small and large s/r Ã . Increasing the strength of selection is inefficient in aiding drift-limited adaptation, in line with the expectation that the effect of genetic drift is only very weakly dependent on selection in two-dimensional habitats [27](see also S1 Fig). This suggests that genetic basis of adaptation is not important for a drift-induced limit to a species' range. Yet it is plausible that there is another limit, in which selection per locus becomes important [27], that arises when the optimum changes abruptly and even when the population (neighbourhood) size is large (i.e., in an entirely different regime). A dedicated synthesis connecting the step-limited and drift-limited regimes would be of a clear interest. Importantly, once genetic drift starts to have an effect, the habitat needs to be fairly broad to be two-dimensional [37]. In narrow habitats (such as in [27]), some dependency of drift-induced expansion threshold on selection per loci would be expected [26]. Note that the apparent independence of the expansion threshold on s/r Ã does not imply that rate of range expansion should also be independent of the strength of selection.
In nature, conditions are unlikely to change uniformly. Abiotic environment (such as temperature, precipitation, solar radiation) does not, in general, change in a linear and concordant manner [38], and neither does the biotic environment, such as the pressure from competitors and predators, which affects the attainable population density and can increase the asymmetry in gene flow [39,40]. I now investigate whether adaptation fails near the expansion threshold as conditions change across space. For example, we can imagine that the population starts well adapted in the central part of the available habitat, and as it expands, conditions become progressively more challenging (see S2A Fig); i.e., the effective environmental gradient B gets steeper. As the expanding population approaches the expansion threshold, adaptive genetic variance progressively decreases below the predicted value [13], which would be maintained by gene flow in the absence of genetic drift (Fig 5A, grey dashed line). This is a result of an increased frequency of demes with limited adaptation, leading to higher rates of extinctions and recolonisations, which reduce both adaptive and neutral diversity (see Fig 5B). Range expansion then ceases at the expansion threshold as the genetic variance drops to the critical value at which only limited adaptation is stable [24], assuming genetic variance is fixed ( Fig  5A, dotted line). This is because although populations can persist with limited adaptation (Fig  4), the transient amount of genetic variance maintained under limited adaptation is almost never consistent with range expansion (see Fig 2, open circles). On a steepening gradient, a sharp and stable range margin forms. This contrasts to uniformly changing conditions (linear gradient, constant carrying capacity) in which populations steadily expand or contract.
In a large population, the ability to adapt to heterogeneous environments is independent of dispersal: this is because both the local genetic variance (measured by standing genetic load), which enables adaptation to spatially variable environments, and the perceived steepness of the environmental gradient (measured by dispersal load) increase at the same rate with gene flow [13]. Yet, in small populations, dispersal is beneficial because the drift-reducing effect of dispersal overpowers its maladaptive effect. This is demonstrated in Fig 6-the neighbourhood size N increases faster with dispersal than the effect of swamping by gene flow (B) does; hence, as dispersal increases, the population gets above the expansion threshold at which uniform adaptation can be sustained. Around the expansion threshold, a small change in dispersal (connectivity) can have an abrupt effect on adaptation across a species' range and the species' persistence. A small increase in dispersal can lead to recovery of uniform adaptation with an This illustrative run shows that as the effective environmental gradient steepens away from the central location, adaptive genetic variance must increase correspondingly in order to maintain uniform adaptation. (a) Median population density stays fairly constant across the range (blue dots), following the deterministic prediction (N, blue dashed line). Genetic variance (black dots) increases due to gene flow across the phenotypic gradient-the deterministic expectation is given by the grey dashed line (see Model section of Methods for details). Yet, as the environmental gradient steepens, genetic variance fails to increase fast enough, and near the expansion threshold, adaptation fails. The dotted line gives the corresponding critical genetic variance, below which only limited adaptation is expected in a phenotypic model with a fixed genetic variance (A % B 2 =2, in which A is the standing genetic load; [24]). (b) As the environmental gradient steepens, the frequency of limited adaptation within the metapopulation increases (black and grey), and hence neutral variation decreases (blue). The black line gives the proportion of demes with limited adaptation after 50,000 generations, when the range margin appears stable; grey gives the proportion after 40,000 generations (depicted is an average over a sliding window of 15 demes). The median is given over the neutral spatial axis Y (with constant optimum); the trait mean, the population trait mean, variance, and population density in two-dimensional space is shown in S3 Fig, which also lists all the parameters. arbitrarily wide continuous range. Further increase of dispersal is merely enhancing the rate of range expansion at the expense of a slight cost to the mean fitness due to rising dispersal load and standing load and can be associated with further costs, such as Allee effect (see, e.g., [17]). Therefore, the expansion threshold provides an interpretation for optimality of an 'intermediate' dispersal, benefiting the species' persistence.

Discussion
Here, I have shown that adaptation fails when positive feedback between genetic drift, maladaptation, and population size reduces adaptive genetic variance to levels that are incompatible with continuous adaptation. The revealed expansion threshold differs qualitatively from the limit to adaptation identified previously [26] for a population living along a one-dimensional habitat. This is because in two dimensions, dispersal mitigates the loss of diversity due to genetic drift more effectively, such that it becomes (almost) independent of selection [34]. The expansion threshold implies that populations with very small neighbourhood sizes (N ≲ 1/2), Open circles indicate limited adaptation, in which a species' range is fragmented and each subpopulation is only matching a single optimum, whilst its genetic variance is very small. As dispersal increases, population characteristics get above the expansion threshold (dashed line), and hence, uniform adaptation becomes stable throughout the species' range. Local population density stays fairly constant, around N = 3.5, whilst total population size increases abruptly above the expansion threshold as the population maintains a wide range (not shown). Parameters for these simulations are given in the Individual-based simulations section of the Methods; the scaling of N and B with dispersal σ is clear from the Methods, section Parameters. The rate of range change is not significantly different from zero for the first three simulations above the expansion threshold; black centre (bottom left) indicates extinction. which suffer a severe reduction in neutral heterozygosity, will be prone to collapse based on demographic stochasticity alone. However, even in the absence of demographic stochasticity, genetic drift reduces the adaptive genetic variance required to sustain adaptation to a heterogeneous environment. The expansion threshold describes when this reduction due to genetic drift is incompatible with continuous adaptation, predicting a collapse of a species' range. If the expansion threshold is reached as the species expands through its habitat, a sharp and stable range margin forms. If there is a drop below the expansion threshold throughout the species' range, as after a sudden drop in carrying capacity, adaptation abruptly collapses throughout a species' range. The result is either extinction or a fragmented metapopulation consisting of a conglomerate of subpopulations, each adapted to a single phenotypic optimum. It follows that near a range margin, we expect increased range fragmentation and a decrease in adaptive genetic variance. The threshold gives a theoretical base to the controversial issue of the importance of evolution (genetics) and ecology (demography) for assessing vulnerability of a species [41,42]. The predicted sharp species' range edge is in agreement with the reported lack of evidence for 'abundant centre' of a species' range, which, although commonly assumed in macroecological theory, has little support in data [3,11,43,44]. Lack of abundant centre is consistent both with uniform adaptation and with limited adaptation in a metapopulation.
The expansion threshold provides a general foundation to species-specific eco-evolutionary models of range dynamics [45]. Its components can be measured in wild populations, allowing us to test the robustness of the theory. First, the effective environmental gradient B can be measured as fitness loss associated with transplant experiments on a local scale, relative to a distance of generational dispersal along an environmental gradient. The environmental gradient can include both biotic and abiotic effects and their interactions [46]-notably, the effective environmental gradient B steepens due to increased asymmetry in gene flow when carrying capacity varies across space, e.g., because of partial overlap with competitors [40]. Second, the neighbourhood size N can be estimated from neutral allele frequencies [47,48]. Estimates of neighbourhood size are fairly robust to the distribution of dispersal distances [49]. Though near the expansion threshold, both the noisiness of the statistics and the homozygosity will increase due to local extinctions and recolonisations [50]. An alternative estimate of neighbourhood size can be also obtained from mark-recapture studies by measuring population density and dispersal (as an approximation for gene flow) independently [47].
Because the expansion threshold is free of any locus-or trait-specific measure, the result appears independent of genetic architecture, readily extending to multiple traits regardless of their correlations (compare to [51][52][53][54][55])-yet the mean fitness will decline because of 'drift load' as the number of traits independently optimised by selection increases [56,57]. Hence, if the fitness landscape is highly complex, the expansion threshold constitutes a lower limit. Naturally, there can be further costs arising in a natural population that I have neglected here, such as the Allee effect [17]. In general, while the numerical constants may change when natural populations deviate in their biology from our model assumptions, the scale-free parameters identified in this study remain core drivers of the intrinsic dynamics within a species' range. Notably, the early classic studies assuming fixed genetic variance [24] predicted that dispersal into peripheral populations is detrimental because it only inflates the effective environmental gradient B. Yet, when genetic variance can evolve, dispersal into small marginal populations also aids adaptation by increasing local genetic variance and by countering genetic drift. The net effect of dispersal into small marginal populations (below the expansion threshold) is then beneficial because their neighbourhood size increases faster with dispersal than the effective environmental gradient B steepens.

Model
I model evolution of a species' range in a two-dimensional habitat, in which both population dynamics and evolution (in many additive loci) are considered jointly. The coupling is via the mean fitness " rð" z;NÞ, which gives the growth rate of the population, and decreases with increasing maladaptation: " rð" z;NÞ ¼ r e ðNÞ þ " r g ð" zÞ. The ecological component of growth rate, r e , can take various forms: here, the regulation is logistic so that fitness declines linearly with density N: r e = r m (1−N/K), in which r m is the maximum per capita growth rate in the limit of the local population density N! 0. The carrying capacity K (for a perfectly adapted phenotype) is assumed uniform across space. The second term, r g ð" zÞ 0, is the reduction in growth rate due to deviation from the optimum. Selection is stabilising: the optimum θ changes smoothly with one spatial dimension (x): for any individual, the drop in fitness due to maladaptation is r g (z) = −(z−θ) 2 /(2V s ). Here, V s gives the width of stabilising selection; strength of stabilising selection is γ = −V P /(2V s ), in which V P = V G +V E is the phenotypic variance. A population with mean phenotype " z has its fitness reduced by " r g ð" zÞ ¼ À ð" z À yÞ 2 =ð2V s Þ À V P =ð2V s Þ. The phenotype z is determined by many di-allelic loci with allelic effects a i ; the model is haploid, hence " z ¼ X i a i p i , in which p i is the allele frequency at locus i. Phenotypic variance is V P = V G +V E . The loss of fitness due to environmental variance V E can be included in r Ã m ¼ r m À V E =ð2V s Þ; V E is a redundant parameter. Selection is 'hard': both the mean fitness (growth rate) and the attainable equilibrium den-sityN ¼ Kr Ã =r m ¼ Kð1 À V G =ð2r m V s ÞÞ decrease with maladaptation. Expected genetic variance maintained by gene flow in the absence of genetic drift is V G ¼ bs ffiffiffiffi ffi V s p [13]; the contribution due to mutation is small, at mutation-section balance V G;m=s ¼ 2mV s n l , in which μ gives the mutation rate per locus and n l the number of loci. Table 1. Three scale-free parameters: B, N, and s/r Ã (top) describe the system. T Middle section gives informative derived parameters. The bottom section gives seven parameters of the model before rescaling, in which the seventh parameter, mutation rate μ, can be neglected because variance maintained by mutation-selection balance, V G , μ/s = 2μV s n l , is typically much smaller than variance generated by gene flow across environments, Is the sky the limit? On the expansion threshold of a species' range

Individual-based simulations
Discrete-time individual-based simulations are set to correspond to the model with continuous time and space. The space is a two-dimensional lattice with spacing between demes of δx = 1. Every generation, each individual mates with a partner drawn from the same deme, with probability proportional to its fitness, to produce a number of offspring drawn from a Poisson distribution with mean of Exp[r(z, N)] (this includes zero). The effective diploid population density N e hence equals half of the haploid population density N, and N = 4πN e σ 2 = 2πNσ 2 . The life cycle is selection ! mutation ! recombination ! birth ! migration. Generations are nonoverlapping, and selfing is allowed at no cost. The genome is haploid with unlinked loci (the probability of recombination between any two loci is 1/2). The allelic effects α i of the loci combine in an additive fashion; the allelic effects are uniform throughout this study, α i α. Mutation is set to μ = 10 −6 , independently of the number of loci. Migration is diffusive with a Gaussian dispersal kernel. The tails of the dispersal kernel need to be truncated: truncation is set to two standard deviations of the dispersal kernel throughout, and dispersal probabilities and variance are adjusted so that the discretised dispersal kernel sums to 1 [58]. Simulations were run at the computer cluster of IST Austria using Mathematica 9 (Wolfram). The code for the simulations, together with a working example, have been deposited as a single Ã .cdf file at Dryad Digital Repository, https://doi.org/10.5061/dryad.5vv37 [36]. This file can be viewed with CDF Player, a free application developed by Wolfram Research, and also contains all the figures with their underlying data.
Parameters. There are in total 10 parameters in the individual-based model, but only 7 are used to describe the model dynamics in continuous time. These are listed at the bottom of Table 1. They are the environmental gradient b = [0.012, 2], dispersal distance σ¼ ½0:1; 1:3, carrying capacity for a well-adapted phenotype K = [3,31], width of stabilising selection V s ¼ ½0:005; 6, the maximum intrinsic rate of increase r m = [0. 2,2], and the mutation rate μ, fixed to μ = 10 −6 . The [min, max] interval gives the parameter range used in the 780 randomly sampled runs, with their distributions described in S3 Fig. The number of genes and demes is not included in the continuous time description (and hence the rescaling) because it assumes that space is not limiting and that all loci have equivalent effect with no statistical associations among them. In the individual-based model, the habitat width is set to be wide enough to be effectively two-dimensional under diffusive dispersal for thousands of generations [37]: 100 dispersal distances σ along the neutral direction and at least 10 cline (deterministic) widths along the gradient. The number of genes contributing to the adaptation across the species' range is n l = [5,2996], with the estimated number of locally polymorphic genes between 1 and 299. Since mutation rate is fixed at μ = 10 −6 , the genomic mutation rate has a wide range, U = [5.10 −6 , 3.

Continuous model
For any given additive genetic variance V G (assuming a Gaussian distribution of breeding values), the change in the trait mean " z over time satisfies: The first term gives the change in the trait mean due to migration with mean displacement of σ; the second term describes the effect of the asymmetric flow from areas of higher density. The third term gives the change due to selection, given by the product of genetic variance and gradient in mean fitness [59,Eq 2]. The last term z gives the fluctuations in the trait variance due to genetic drift: z ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi V G;LE =N q dW z ðx;y;tÞ, in which dWÃ represents white noise in space and time [34,60]. V G;LE ¼ X i a 2 i p i q i denotes genetic variance assuming linkage equilibrium. The trait mean is " z ¼ X i a i p i for a haploid model, in which p i is the i-th allele frequency, q i = 1−p i and α i is the effect of the allele on the trait-the change of the trait mean " z as frequency of locus i changes from 0 to 1. For both haploid and diploid models, the allele frequencies p i change as: The expected change of allele frequency due to a gradient in fitness and local heterozygosity is p i q i @" r @p i ¼ s i p i q i ðp i À q i À 2D i Þ, in which selection at locus i is s i a 2 i =ð2V s Þ and D i ¼ ð" z À bxÞ=a i [13,Appendix 3]. Here, the fourth term describes the change due to (symmetric) mutation at rate μ. The last term ε describes genetic drift [34, Eq 7]: ε ¼ ffiffiffiffiffi p i q i N q dW ε ðx;y;tÞ, in which N is the haploid population density. Population dynamics reflect diffusive migration in a two-dimensional habitat, growth due to the mean Malthusian fitness " r, and stochastic fluctuations. The number of offspring follows a Poisson distribution with mean and variance of N; fluctuations in population numbers are given by [61]: x ¼ ffiffiffiffi N p dW x ðx;y;tÞ:

Continuous model: Rescaling
The model can be simplified by rescaling [13,59] time t relative to the strength of density dependence r Ã , distance x relative to dispersal σ, trait z relative to strength of stabilising selection 1=ð2V s Þ and local population size N relative to equilibrium population size with perfect adaptation:N ¼ Kr Ã =r m , T = r Ã t; X ¼ x ffiffiffiffi ffi 2r Ã s 2 q ; Z ¼ z ffiffiffiffiffiffi r Ã V s p ;Ñ ¼ N=N . Note that near the equilibrium of a well-adapted population,Ñ % 1.
environmental gradient steepens, so does the predicted genetic variance maintained by gene flow in the absence of genetic drift (grey mesh; V G ¼ bs ffiffiffiffi ffi V s p , [13]). Past the threshold, genetic variance starts to abruptly fall off the prediction, and adaptation begins to fail. The median values (across the neutral Y-space) are shown in , with median s/r Ã = 0.006. Not pictured is the uniform mutation rate per locus (μ = 10 −6 ), the number of demes along the two spatial directions, which was set to be at least 10 deterministic cline widths wide along the spatial axis with environmental gradient (X), whereas along the neutral spatial axis (Y), the habitat width was 100σ. (EPS)