Skip to main content
Open Access

Local Adaptation Interacts with Expansion Load during Range Expansion: Maladaptation Reduces Expansion Load


The biotic and abiotic factors that facilitate or hinder species range expansions are many and complex. We examine the impact of two genetic processes and their interaction on fitness at expanding range edges: local maladaptation resulting from the presence of an environmental gradient and expansion load resulting from increased genetic drift at the range edge. Results from spatially explicit simulations indicate that the presence of an environmental gradient during range expansion reduces expansion load; conversely, increasing expansion load allows only locally adapted populations to persist at the range edge. Increased maladaptation reduces the speed of range expansion, resulting in less genetic drift at the expanding front and more immigration from the range center, therefore reducing expansion load at the range edge. These results may have ramifications for species being forced to shift their ranges because of climate change or other anthropogenic changes. If rapidly changing climate leads to faster expansion as populations track their shifting climatic optima, populations may suffer increased expansion load beyond previous expectations.

Online enhancements:   appendixes, video files. Dryad data:


Species range expansion is a complex process that is widely studied in evolutionary biology, ecology, and conservation biology (Hastings et al. 2005; Phillips et al. 2006; Excoffier et al. 2009; Hallatschek and Nelson 2010; Chen et al. 2011; Colautti and Barrett 2013). An array of ecological and evolutionary factors are known to affect the success or failure of range expansion, such as dispersal limitation (Hastings et al. 2005; Marsico and Hellmann 2009; Hargreaves et al. 2014), interspecific competition (Case and Taper 2000; Price and Kirkpatrick 2009; Svenning et al. 2014; Louthan et al. 2015), or an inability to adapt to new conditions (Angert et al. 2008; Holt and Barfield 2011; Polechová and Barton 2015). A recently growing field of research has focused on an additional effect: the genetic load accumulated during the expansion process (Excoffier et al. 2009; Hallatschek and Nelson 2010; Peischl et al. 2013, 2015; Peischl and Excoffier 2015). In particular, the increasing frequency of deleterious variants in human populations that have expanded out of Africa has been widely studied (Lohmueller et al. 2008; Simons et al. 2014; Do et al. 2015; Henn et al. 2015b). A further process that can increase genetic load during range expansion is maladaptation to the local environment due to gene flow from foreign environments, known as migration load (Kirkpatrick and Barton 1997; Barton 2001; Polechová and Barton 2015). In this study, we explore the impacts and interactions of heterogeneous selection due to an environmental gradient and the accumulation of deleterious mutations during range expansion.

During range expansions, an intriguing suite of population genetic processes can change the course of evolution compared to that in populations that are growing in size without movement. Populations at the expanding front of a species range undergo serial founder events, where each new colonization into further territory creates a population bottleneck, leading to reduced genetic diversity in what becomes the new edge population and the primary source of future colonists. These processes create a persistently reduced effective population size at the range edge, decreasing the efficacy of selection and increasing the strength of random genetic drift. This leads to the process of allele surfing (Klopfstein et al. 2006), whereby a deleterious mutation that arises at the range edge is more likely to drift to high frequency than it would if it arose in the denser core of the species range. Because selection is inefficient, strongly deleterious alleles may persist and reach higher frequencies than expected in a population at equilibrium (Peischl and Excoffier 2015). Theoretical work shows that allele surfing is the main cause of increased deleterious allele frequency within recently expanded populations (Excoffier et al. 2009). The reduction in fitness due to the accumulation of these deleterious alleles is termed expansion load (Peischl et al. 2013; Peischl and Excoffier 2015).

There is controversy over empirical evidence for the accumulation of deleterious mutations of both higher frequency and increased effect sizes in recently expanded populations. Expansion load has been posited as an explanation for deleterious alleles in human populations that have undergone recent geographic expansions, such as genetic diseases in populations of Quebec (Labuda et al. 1997; Scriver 2001; Yotova et al. 2005) and Scandinavia (Norio 2003). Whether the effects of demographic history and expansion load are great enough to accumulate significant levels of deleterious mutations to affect population fitness and persist into future generations of humans is still debated (Lohmueller 2014a, 2014b; Simons et al. 2014; Sousa et al. 2014; Gravel 2016). Yet it has been shown that, as human expansion out of Africa proceeded, deleterious mutations of larger effect rose in frequency, shifting the entire distribution of allelic effect sizes upward through time (Henn et al. 2015b). Expansion load may thus have serious repercussions in many other species that undergo expansions, such as those recolonizing formerly glaciated areas (Hewitt 1999), colonizing a new continent (Sakai et al. 2001), or tracking recent climate change (Chen et al. 2011). Understanding the detriment caused by this process may aid in efforts to combat invasive species or in assisted migration. As global climate change proceeds, more species may be required to move their ranges in order to survive, and if fitness losses simply due to the expansion process are common, populations may be left at higher risk for extinction from stochastic, catastrophic events.

A second evolutionary process that can affect range expansion is local adaptation to heterogeneous environmental conditions. Selective environments can vary over space: for example, environmental gradients in temperature or photoperiod from equatorial to polar latitudes influence important traits (Conover 1992; Montague et al. 2008). A foundational study by Kirkpatrick and Barton (1997) showed that steep environmental gradients can impede range expansion because of the migration of individuals from the range core to the range edge, preventing local adaptation to edge conditions. Central populations existing on an environmental gradient receive symmetric migration from populations both higher and lower on the gradient, leaving the population mean unchanged. In contrast, at a range edge, migration is asymmetric from core to edge, as small edge populations receive proportionally more immigration from the denser species core, causing edge populations to be swamped by locally maladaptive alleles and diverge from their local optimum. Therefore, when migration rates are high enough, the edge population can fail to adapt. When the environmental gradient is steep enough, the edge populations can experience sufficient fitness reductions to result in local extinction and prevent range expansion. Further studies (Barton 2001; Polechová and Barton 2015) have increased the biological realism of range expansion models (including evolution of genetic variance and effects of genetic drift), confirming theoretically that evolutionary processes alone can lead to the formation of a stable range edge.

It is clear that expansion load can result in reduced fitness of recently expanded populations and that heterogeneous selection along an environmental gradient can result in local maladaptation of populations at the edge of a species range. These two evolutionary processes—maladaptation from underlying environmental gradients and expansion load from surfing of deleterious mutations at range edges—may occur simultaneously in many range expansions and could interact in interesting ways that have not previously been studied. We consider two contrasting a priori hypotheses for how they may interact. One possibility is that local maladaptation and expansion load may interact positively (i.e., each making the effects of the other more pronounced). This may be expected because both processes could lead to smaller population sizes at the expanding front, and both are stronger with greater genetic drift: expansion load increases because greater drift would cause more deleterious alleles to increase in frequency, and maladaptation worsens because smaller populations receive disproportionately more immigration and greater drift would eliminate more of the genetic variance needed to respond to local selection. Alternatively, expansion load and local maladaptation could interact negatively, with each partially ameliorating the effects of the other. This is plausible because both forms of load could reduce the pace at which ranges expand by reducing the fitness of individuals at the range margin. Slower expansion enables fitter alleles to reach the edge through migration and to rescue populations suffering from accumulated expansion load. Similarly, slower range expansion could allow more time for populations to adapt to local conditions.

We investigate these alternative hypotheses to establish the role of expansion load and migration load in isolation as well as in combination. Our goal is to employ the most biologically reasonable parameters possible within the realm of our model. Using individual-based simulations on a two-dimensional, spatially explicit, and approximately continuous landscape, we compare the reductions in fitness (load) that populations experience during range expansions over a series of environmental gradients and deleterious mutation rates. Our results have implications for the predicted prevalence of expansion load in species across various types of environments and inform our understanding of the complex demographic and genetic processes that occur during range expansions.


We model a species range expansion forward in time, using the simulation program Nemo (Guillaume and Rougemont 2006). The population undergoes an initial burn-in period to mutation-selection equilibrium, after which it is allowed to expand across an empty landscape. Generations are nonoverlapping. Individuals are monoecious and diploid, possessing both a quantitative trait experiencing stabilizing selection and a set of loci subject to unconditionally deleterious mutations. We compare three environmental gradients over which expansion occurs and three genome-wide deleterious mutation rates. All combinations of parameter values are listed in table A1. We ran 20 independent replicates for each of our simulated scenarios, with the data available from the Dryad Digital Repository: (Gilbert et al. 2017).

The Model

We model space explicitly on a 2,000×40-unit landscape, using a modified version of Nemo available at Each 1×40-unit slice of the landscape is referred to as a cross section (fig. 1). This modification allows a discrete approximation of continuous space. The landscape is divided into 80,000 square units, which we call cells. Within a cell, space is undefined and all individuals are considered to exist simultaneously at the center of the cell. The dispersal kernel and breeding window occur across cells, allowing individuals to interact over a larger spatial scale. Life-cycle events occur in the order breeding, dispersal, selection, and then population regulation. Subpopulation regulation enforces a carrying capacity of 6 individuals per cell, above which individuals are randomly culled at the end of each generation. Since breeding is not limited to within one cell, the realized neighborhood size, as defined by Wright (1946), is approximately 300 individuals.

Figure 1.
Figure 1.

Small-scale representation of the landscape simulated with a grid of cells where the burn-in period occurs in the landscape core at the left, expansion proceeds to the right, and any vertical column of cells is a cross section. Actual landscapes were 2,000×40 cells, with a 40×40-cell core.

Populations were initiated at carrying capacity in the leftmost 40 cross sections of the landscape, which we term the landscape core (fig. 1). After a burn-in period of 15,000 generations, the remaining 1,960 cross sections of the landscape become available, allowing for expansion to occur. The environmental optimum is constant across any given cross section.

The environmental gradient underlying the landscape changes in only one dimension, along the axis of expansion. We compare three values for the steepness of this gradient, b, which defines the change in phenotypic optimum over space: a flat landscape with no gradient (b=0.0), a shallow gradient (b=0.0375), and a steep gradient (b=0.375). On the shallow gradient, moving the average dispersal distance along the X-axis (1.596) would result in a fitness loss of 0.02% if an individual were perfectly adapted in its natal environment, whereas on the steep gradient, an average dispersal event would reduce an individual’s fitness by 2.4%, given perfect adaptation in its natal environment. Each of these gradients results in a different equilibrium level of migration load in populations.


Our modified version of Nemo uses a breeding window within which individuals may search for a mate within their own cell on the landscape or in nearby cells. For a single mating event, each individual searches for a mate on the basis of its relative mating probability. The absolute mating probability follows an approximate bivariate Gaussian function where f(Δx,Δy)exp[(Δx2/2σbreed2+Δy2/2σbreed2)] gives the relative probability of choosing a mate within a cell at distance Δx and Δy along the axes of the landscape away from the natal cell; σbreed defines the size of this breeding window and was held constant at 0.5 landscape units, and the maximum searchable distance was restricted to 4σbreed. Because space within a cell is undefined, we arbitrarily assign a cell width and integrate the mating probability over both dimensions of each cell within the breeding window, based on the distance to the center of the natal cell, to discretize the continuous probability distribution to the relative probability of selecting each individual in the given cell with which to mate. The probability for each cell is then multiplied by the number of potential mates in the cell, and these products are summed over all cells. The absolute probability of selecting each individual is then that individual’s probability divided by this sum over all individuals within the range of the breeding window. This results in a two-dimensional breeding window in which a potential mate in the mating individual’s focal cell is the most likely to be chosen and those in more distant cells (of the 12 surrounding cells) are chosen with decreasing probability. An individual self-fertilizes only when no other individuals are present within the breeding window. This will most often happen at the expanding range edge, where population densities are lowest, and mimics the ability of many plant species to self under conditions of pollen limitation (Hargreaves and Eckert 2014). Each female’s fecundity was drawn from a Poisson distribution with mean 7 to determine the number of mating events.


We modeled dispersal for most simulations according to an approximate bivariate Gaussian kernel, similar to the breeding kernel: f(Δx,Δy)exp[(Δx2/2σdisperse2+Δy2/2σdisperse2)]. The dispersal kernel gives the forward (in time) migration probabilities for offspring to disperse to a given patch distance Δx and Δy away. The maximum dispersal distance was capped at 8σdisperse units. The “standard deviation” of the Gaussian dispersal kernel σdisperse was set to equal 2 landscape units. We also simulated a leptokurtic kernel to test the effects of rare, long-distance dispersal events. The leptokurtic dispersal kernel was a mixture distribution created from a weighted sum of two different Gaussian kernels (Ibrahim et al. 1996). This summed a first kernel having σdisperse=1.5 with a second kernel having σdisperse=5.5, where the first kernel was weighted by 0.89 and the second by 0.11. This created our leptokurtic kernel with the same average dispersal distance as the Gaussian kernel and a total value of kurtosis equal to 10, which is not unrealistic for long-distance dispersal in many species of plants and animals (Skalski and Gilliam 2000; Lowe 2009; Guttal et al. 2011). The leptokurtic kernel was truncated to have the same maximum dispersal distance, which in both cases meant that an individual could disperse at most 16 units in any one direction from its natal cell. Borders of the 40×2,000-unit landscape were absorbing, so any individual migrating beyond the edge was removed. R code for calculating and discretizing the dispersal and breeding kernels is available in a wrapper package written for Nemo, aNEMOne, available at


We modeled a genetic architecture where 100 quantitative trait loci and 1,000 loci subject to unconditionally deleterious mutations were randomly and independently placed on the genome for each simulation. Each genome consisted of 10 chromosomes of 100 cM each, where loci separated by 1 cM have a 1% chance of a crossover event between them each generation. We selected realistic mutational parameters for these loci as follows.

The quantitative trait, z, was controlled by 100 additive loci. Each allele at these loci can take any real value. In this model an individual’s trait value is given by the sum of allelic effects across all quantitative trait loci, with no dominance. Alleles are continuous, such that a mutation’s effect is added to the existing allelic value. The component of fitness (survivorship) for the quantitative trait value z is wz=exp[(zzopt)2/2ω2], where zopt is the optimal trait value for the given location on the landscape. At the start of the burn-in, all individuals were initiated with z equal to the mean zopt of the burn-in area (core).

The evolution of trait z under stabilizing selection depends on the mutational variance, VM, and the inverse strength of stabilizing selection on the trait, ω2. These properties have been empirically estimated for a number of traits, along with the corresponding environmental variances, VE. For our model we set VM and ω2 relative to the same arbitrary VE value (VE=1). There is evidence that ω25VP, on average (where VP is the phenotypic trait variance; Kingsolver et al. 2001; Johnson and Barton 2005). The relationship between VP and VE can be expressed in terms of heritability, h2=1VE/VP, and gives a value of ω25VE/(1h2). Given VE=1 and a typical heritability of h21/3 (Mousseau and Roff 1987; Houle 1992), we set ω2=7.5.

Most reported values of VM/VE are in the range of 10−4 to 10−3 (Houle et al. 1996). The value of VM depends on the genome-wide rate of mutations affecting the trait, Uz, and the expected squared effect of a mutation on the trait, E[α2], where VM=UzE[α2]. If mutational effects are Gaussian with E[α]=0, then E[α2]=V[α], the variance of mutational effects on the trait. These components are difficult to estimate, but according to one theoretical approach we expect VPVE=4Uzω2 (Turelli 1984; Charlesworth and Charlesworth 2010). Given the parameters above, this implies Uz0.02, similar to some direct estimates (Lynch and Walsh 1998). In our simulations we used 100 quantitative trait loci, with a mutation rate per locus of 10−4, giving a diploid mutation rate, Uz, of 0.02. We further assumed V[α]=0.02 (i.e., mutational effects on z are drawn from a normal distribution with mean 0 and variance 0.02), giving a value for VM of 4×10−4.

We also modeled 1,000 biallelic loci subject to unconditionally deleterious mutations. Fitness for these mutations was multiplicative across loci, given by wD=i=11,000(1hisiψhet,isiψhom,i), where si and hi are the selection and dominance coefficients, respectively, for the deleterious allele at locus i and ψi indicates the presence (1) or absence (0) of a deleterious allele at locus i in the heterozygous or homozygous state, respectively. An individual’s fitness is the product of wz and wD.

We considered genome-wide diploid deleterious mutation rates, UD, of 0.1 and 1.0 in order to encompass probable rates for a variety of taxa, including Drosophila melanogaster (Haag-Liautard et al. 2007), Caenorhabditis elegans (Denver et al. 2004), Arabidopsis thaliana (Shaw et al. 2000), Amsinckia sp. (Schoen 2005), and possibly nonhuman endothermic vertebrates (Baer et al. 2007). However, we note that in humans UD likely exceeds 1.0 (Keightley 2012), and UD is likely less than 0.1 in many other organisms (Baer et al. 2007; Halligan and Keightley 2009). We also included a treatment where deleterious mutations were absent (UD=0.0). Deleterious-mutation rates per haploid locus were thus 0, 5×10−5, and 5×10−4, and we allowed for back-mutation at rates of 0, 5×10−7, and 5×10−6, respectively. Once mutated, a locus could not further mutate, and the only possible change would be due to a back-mutation. These 1,000 loci do not accurately portray the number of possible loci in biological systems but are instead an approximation required by the constraints of simulation. Because we matched our genome-wide mutation rates to realistic values, each of these loci is more representative of a region of the genome within which a deleterious mutation may arise. Thus, for the distribution of effects we use, described below, small-effect mutations may be approaching saturation, but because these small-effect loci do not contribute substantially to fitness reductions, this lack of realism in terms of the number of loci is not detrimental to the accuracy of our model.

Although the true distribution of mutational fitness effects is empirically difficult to measure and may be complex, there is evidence that most mutations have small effects on fitness and that rare mutations have large effects on fitness (Eyre-Walker and Keightley 2007). We modeled the homozygous fitness effects (s) of deleterious mutations by using a leptokurtic gamma distribution with mean 0.01 and shape parameter 0.3, such that most mutations have s<1% (Keightley 1994). The mutational effect of each locus was drawn from this distribution at the start of each independent simulation run and remained constant throughout the run.

Estimates of the average dominance (h) of deleterious mutations are scarce, range widely, and are subject to various biases (Halligan and Keightley 2009; Agrawal and Whitlock 2011). It is likely that most new mutations are partially recessive, and there is evidence for a negative relationship between h and s (Agrawal and Whitlock 2011). We follow previous authors (Lynch et al. 1995; Deng and Lynch 1996) in assuming an exponential relationship between h and s, as well as a mean h of approximately 0.37. Specifically, we assume h=exp(−51.1s)/2. In addition, we included a class of deleterious mutations with s=1 and h=0, that is, recessive lethals, and assumed that 3% of deleterious mutations fall into this class, to match the genome-wide rate typically observed in D. melanogaster (Fry et al. 1999).

Mimicking Mutation Load

To disentangle the effects of mutation load and expansion load, we compared an additional parameter set to the core set described above. When populations are at mutation-selection balance, the reduction in fitness due to deleterious mutations is termed the mutation load. During range expansion, deleterious mutations can increase because of genetic drift at the expanding front. This excess load beyond mutation load is termed the expansion load. To mimic the presence of mutation load only, without expansion load, we set UD to 0 and altered individual fecundity to reduce realized fitness to a level that matches the mutation load measured in the range core, which has not undergone expansion, for each of our scenarios UD=0.1 and UD=1.0. This was done only for the cases of Gaussian dispersal. We calculated the equilibrium fitness due only to deleterious mutations in the core landscape populations and reduced mean fecundity to recreate this same realized fitness. Mean fecundity in the case matched to UD=0.1 was approximately 6.475, and that for the UD=1.0 case was approximately 3.721.


We assess the impact of expansion across an environmental gradient and in the presence of deleterious mutations both independently and in combination. We quantify two measures from the simulation results: the speed of range expansion and mean fitness at the range edge versus that in the core. Fitness is measured after population regulation and therefore includes only individuals surviving each generation. We partition fitness into each of its contributing components: wz for the quantitative trait and wD for the deleterious alleles.

To track the expanding front, we use a single landscape cross section as the range edge. This edge was defined as the first cross section away from the core at which population size is at 50% of the core’s equilibrium population size. The speed of expansion measures the number of cross sections over which this edge cross section travels per generation from the end of the burn-in until reaching the second-to-last landscape cross section.

We examine fitness at the range edge by using a broader definition of the edge to reduce sampling error. We defined this measure of the range edge to contain all individuals present within the cross section used to measure expansion speed as well as all individuals present in cross sections farther toward the empty landscape. Expansion load measured the excess load that accumulated at the range edge beyond mutation load, calculated as


Because expansion proceeded at different rates in different scenarios, to create an equivalent comparison of fitness on the landscape we measured populations at the latest recorded point in the simulations before any individuals had dispersed into the last 20 cross sections of the landscape. Since individuals can disperse at most 16 units, this prevented bias from edge effects that might arise upon filling the landscape. We averaged population sizes and fitness within cross sections of the landscape, as we saw no major variation along the axis perpendicular to expansion.

We also examined the distribution of effect sizes for deleterious mutations accumulated at the range edge. We examined the distribution of alleles contributing to expansion load across cases of UD and b. Deleterious mutations were binned by effect size into 40 quantiles based on the underlying gamma distribution of homozygous effects (resulting in approximately 25 loci per bin). Within a given simulation, we calculated the average load due to loci present within each of these bins as follows. For the edge or the core, we calculated fitness in a given bin of n loci as wD=i=1n(1hisiɸhet,isiɸhom,i), where ɸhet, i is the frequency of heterozygotes for a deleterious allele at a given locus i and ɸhom, i is the frequency of homozygotes. Expansion load was then calculated for each bin with equation (1). Because fitness is multiplicative across loci, we then calculated average expansion load per locus in each bin as

.11expansion loadn.


Expansion Speed

No parameter combinations that we investigated ever led to the formation of a stable range edge. However, the rate at which populations expanded across the landscape was affected by both heterogeneous selection from the environmental gradient and deleterious mutations. Steeper environmental gradients slowed expansion (fig. 2; table A1). Compared to the case of no gradient, the steep gradient slowed expansion in the Gaussian-dispersal cases by 77%–80%. Leptokurtic dispersal kernels resulted in faster expansion than Gaussian kernels over the landscape, particularly in the absence of an environmental gradient, where expansion speed was 32%–52% greater for leptokurtic dispersal. Increasing the deleterious genomic mutation rate decreased the speed of expansion. Compared to UD=0, UD=1 for the Gaussian kernel decreased expansion speed by 31.6% with no gradient and by 21.3% on the steep gradient. Cases without expansion load (fecundity adjusted) expanded faster than cases with expansion load at both UD values only in the absence of a gradient, while all cases with a gradient exhibited slight increases or decreases in speed. The greatest speed increase from removing expansion load was by 18.5% for UD=1 and b=0.

Figure 2.
Figure 2.

Average speed of expansion across scenarios. Circles indicate Gaussian dispersal kernels, while triangles indicate leptokurtic dispersal kernels. Dashed lines indicate fecundity-adjusted simulations with Gaussian dispersal, which replicate the loss in realized fitness due to mutation load but lack expansion load because no deleterious alleles are present. The 95% confidence intervals are smaller than the plot points in all cases, and points are offset for better visualization.

Fitness and Load

In general, heterogeneous selection and deleterious mutations both reduced mean fitness in the core and edge, as expected. Furthermore, fitness at the expanding edge was always reduced relative to that in the core of the species range (fig. 3). Overall fitness at the range edge correlated with expansion speed, but the presence of an environmental gradient more strongly changed the speed of expansion (fig. A1). The quantitative fitness component, w¯z, in the core was nearly identical across UD cases and ranged from 0.87 to 0.93 (table A1). For the quantitative trait, fitness reduction at the edge occurred only in the presence of an environmental gradient and ranged from a 25%–77% reduction in w¯z at the edge relative to the core for the Gaussian-dispersal cases. The most severe fitness loss at a range edge (77%) was for UD=0 on the steep gradient (b=0.375).

Figure 3.
Figure 3.

Average core and edge fitness components for the quantitative trait and deleterious alleles. Dashed lines indicate fecundity-adjusted simulations, which replicate the loss in realized fitness due to mutation load but lack expansion load because no deleterious alleles are present. Error bars indicate 95% confidence intervals, and points are offset for better visualization.

The fitness component for loci with deleterious alleles, w¯D, in the range core reflected the equilibrium mutation load reached in the respective UD cases. For UD=0.1, core w¯D was 0.92, while for UD=1, core w¯D was 0.53. Edge w¯D ranged from 0.80 to 0.89 for UD=0.1 and from 0.32 to 0.44 for UD=1. Table A1 reports all fitness values for the respective scenarios simulated.

Expansion load for our parameter sets caused at most a 39% reduction in fitness, while at its weakest it resulted in a 3.5% decrease in fitness (fig. 4). The UD=1 case had a 2.95-fold increase in expansion load over the UD=0.1 case in the absence of an environmental gradient, while on the steepest gradient this change in UD resulted in a 5.1-fold increase in expansion load.

Figure 4.
Figure 4.

Average expansion load across all combinations of environmental gradients and genome-wide deleterious mutation rates. Error bars indicate 95% confidence intervals.

Leptokurtic dispersal increased expansion speed but did not significantly change fitness from the Gaussian results (fig. A2). At a given point on the landscape, fitness recovered after the expanding front passed (fig. A3). Supplemental movies (videos B1B3) show the effects on fitness reduction at the range edge due to surfing and local maladaptation as well as recovery through time due to immigration from the range core.

Interaction of Expansion Load and Heterogeneous Selection

We next investigate how each individual component of fitness is affected by reduced fitness in the other component. We examine the impact of expansion load on the level of local maladaptation in populations and, alternatively, whether local maladaptation affects the degree of expansion load. First, we find that increased load due to deleterious alleles improves the level of adaptation to the local environment at the range edge (fig. 3a). This becomes clear when we consider the fecundity-adjustment simulations that lack expansion load. Within a given value of the environmental gradient, w¯z increases at the expanding front when mutation load is the only effect present (fecundity-adjusted runs) and increases even further in the presence of both mutation load and expansion load (UD=0.1 and UD=1). The largest increase in quantitative trait fitness occurs on the weak gradient (b=0.0375) between UD=0 and UD=1, where the presence of mutation load leads to a 30% improvement in w¯z on average and the presence of mutation and expansion load together increases w¯z by 70%.

Second, we find that the severity of expansion load is substantially reduced by the presence of an environmental gradient (fig. 4). Expansion load decreases with increasing gradient steepness. In the case of UD=1, there is a 54% reduction in expansion load from no gradient to the steepest gradient. In the case of UD=0.1, expansion load is reduced by 73% on the steepest gradient, compared to the case with no environmental gradient.

Loci Contributing to Expansion Load

We investigated the distribution of effect sizes for deleterious mutations accumulating on the expanding population front. For UD=1, the majority of expansion load is due to alleles of intermediate effect (fig. 5b). On the steep gradient (b=0.375), there is less load overall due to each locus but a more even distribution across loci of both intermediate and relatively large effect (e.g., s>0.04). On the shallow gradient and with no gradient, there are fewer large-effect alleles contributing to load relative to the steep gradient but also many more intermediate-effect alleles contributing to more overall load. For UD=0.1 (fig. 5a), there is less overall load and so lower averages within most ranges of allelic effect size. However, we see an inflation of larger-effect alleles above that seen for the UD=1 case in the absence of an environmental gradient.

Figure 5.
Figure 5.

Average expansion load per locus at the range edge. Loci are binned into quantiles drawn from our modeled gamma distribution, shown by the histogram of loci density. The leftmost bins exceed the height of the Y-axis. Average load per locus is then calculated for each of these bins and plotted for each value of the environmental gradient.

There are very few loci present in the larger-effect-size classes, and for UD=0.1 and b=0 these loci contribute disproportionately more to expansion load (fig. 5a). As can be seen in figure 6, some of these larger-effect loci have fixed at the range edge (see also fig. A4). The effect of surfing that contributes to the fixation of these strongly deleterious alleles can be seen in animation over time in video B4, where recovery from fixation follows behind the expanding wave front.

Figure 6.
Figure 6.

Deleterious allele frequencies across the landscape, by homozygous effect size, s. This example is at 250 generations after burn-in for the case of a leptokurtic dispersal kernel, with UD=0.1 and b=0. This example was not randomly chosen but was selected to show that large-effect loci can locally fix under these conditions. The locus indicated by the thickest black line has an effect size of 0.0898.


Range expansions are unique demographic events that lead to an interesting suite of population genetic processes. These processes have been widely studied, yet the combination of both a heterogeneous environmental gradient and expansion load from deleterious mutations has not previously been investigated. Previous theoretical studies focusing independently on either the expansion load or adaptation along environmental gradients predicted reduced fitness at the range edge. Our main result is that expansion load and the load due to local maladaptation interact, so that the load at the expanding range edge is not as great as would be expected from a simple combination of the two. Local maladaptation at the range edge is not as severe in the presence of expansion load, and, of an even greater effect, expansion load is not as severe in the presence of a strong environmental gradient causing local maladaptation. Whether fitness is reduced because of local maladaptation or because of expansion load slows the rate of range expansion by different amounts. This affects both the degree of genetic drift occurring at the range edge and the amount of migration reaching populations at or near the expanding front contributing to the magnitude of load that accumulates during expansion. This interaction is important for predicting the dynamics of modern range expansions due to natural phenomena or human-induced climate change.

We believe that the dominant reasons differ for these two patterns of interaction. First, let us consider the improvement in the local adaptation at the range margin in the presence of expansion load. For a population to expand, the individuals at the expanding front must have an absolute fitness greater than 1, on average. This means that a range margin will occur near where absolute fitness drops below 1, as such populations are sinks persisting by immigration. As a consequence, when there is greater load in one component of fitness, such as that caused by deleterious mutations across the genome, the population will not persist unless other fitness components are great enough to allow sufficient absolute fitness. Therefore, when expansion load is greatest, as occurs at the range edge, the fitness due to local adaptation must be higher in order for the population to have a large-enough fitness to even exist. Note that this improvement in local adaptation is more a case of expansion load eliminating maladapted populations at the expanding front than a reduction in load. As a result, greater adaptation in the quantitative trait is found in the presence of expansion load.

Second, we also observe that expansion load can be greatly reduced in the presence of local maladaptation. In this case, we believe that there is an additional factor explaining this interaction. In the absence of an environmental gradient, where the greatest expansion load accumulated (39%), range expansion is mainly limited by dispersal ability. However, with local maladaptation caused by an environmental gradient, expansion is further slowed by the need for colonizing populations to adapt to the novel local environment. In fact, range expansion is slowed considerably more by a changing selective environment than by a high frequency of unconditionally deleterious alleles (see fig. A1). As a result, the local maladaptation caused by a heterogeneous environment causes the rate of range expansion to slow substantially (fig. 2), which allows more time for migrants to reach marginal populations and less time for deleterious mutations to increase as a result of drift. As range expansion slows, selection has more opportunity to reduce the frequency of deleterious alleles at the edge. Moreover, with slower expansion, high-fitness alleles from the center of the range have more time to disperse toward the edge, restoring genetic diversity at sites locally fixed for deleterious alleles by previous drift, and back-mutation also has more time to generate beneficial diversity.

The mechanism altering expansion speed is a vital part of the process for these interactions to occur. When expansion speed is changed as a result of altered dispersal abilities, there is no expectation for faster expansion to contribute to a further increased load. This explains the lack of further fitness reductions from increased expansion speed in our simulations with long-distance dispersal. Instead, local adaptation at the edge was slightly reduced and expansion load less severe relative to the Gaussian-dispersal case, suggesting that this difference in dispersal models led to only slightly different amounts of expansion load accumulating, given the different amount of connectivity between the core and the edge. Furthermore, we do not interpret the differences in expansion speed as reflecting any innate differences in species range sizes, as we expect range size to instead relate to species characteristics such as dispersal ability and adaptation over many traits to different environmental aspects that we may not have modeled.

The effect sizes of mutations underlying expansion load have important empirical implications, as many studies today aim to understand expansion load in humans after expansion out of Africa (Lohmueller 2014a, 2014b; Henn et al. 2015a, 2015b; Gravel 2016). Whether expansion load exists but is difficult to detect has been debated, and a full understanding of the selection coefficients across human (or other) genomes is lacking (Hancock et al. 2011; Lohmueller 2014b; Simons et al. 2014; Henn et al. 2015a, 2015b). We find an interesting difference in the makeup of expansion load between our two simulated cases of genome-wide deleterious mutation rates. While deleterious alleles of moderate and large effect are relatively rare, we find that they are responsible for the majority of expansion load. These alleles, which would otherwise tend to be purged or kept at low frequency by purifying selection, are able to increase in frequency on the expanding wave front. Interestingly, we find that in cases exhibiting faster range expansion, a larger proportion of expansion load comes from alleles of moderate and large effect, because with faster expansion there will be greater genetic drift at the margin, allowing the rise in frequency of strongly deleterious alleles.

Implications, Caveats, and Future Directions

Our finding that local maladaptation interacts with expansion load has broad evolutionary and ecological implications, including for studies of natural range expansions under climate change, invasive species, and conservation efforts (Hunter and Hutchinson 1994). For example, highly invasive species are known for their rapid rate of spread (e.g., cane toads in Australia; Phillips et al. 2006), providing interesting opportunities to examine whether and how these species accumulate expansion load. Conservation efforts that aim to reintroduce genetic diversity through assisted migration would clearly benefit edge populations in terms of reducing expansion load, but potential disruption to local adaptation is also a concern (Aitken and Whitlock 2013). Finally, both expansion load and migration load are likely to affect climate-change-induced range shifts. Rapid climate change would necessitate fast range expansion. Therefore, expansion load may be increased and local adaptation decreased, leaving struggling populations subject to stochastic extinction events. If the speed of climate change is not too fast, however, populations adapting as they move over space may reduce any potential impacts of expansion load.

Even though our simulation design differs from previous models investigating expansion load, we still find substantial load accumulating throughout the course of range expansion. The presence of hard versus soft selection can change the amount of expansion load accumulated (Peischl et al. 2013; Peischl and Excoffier 2015), as can one- versus two-dimensional landscape models (Peischl et al. 2013). Several other differences, such as a range of mutational effects for deleterious alleles and continuous dispersal rather than a stepping-stone model, might also contribute to the amount of expansion load.

There are several biological features of organisms and their environments that we did not consider in our simulations because of the vast computational resources already required, but these merit future investigation. Our simulations allowed for individuals to self-fertilize when mates were limited, but this is not possible in many species. The inability to self-fertilize could slow range expansion, leading to a reduction in expansion load. Other effects can similarly reduce fitness in small edge populations. Allee effects (Taylor and Hastings 2005) or aggregating dispersal behavior that discourages colonization of empty habitat (Altwegg et al. 2013) would slow expansion, as might dispersal barriers or increased encounters with antagonistic species (e.g., competitors, pathogens; Case et al. 2005; Kubisch et al. 2013). It would be interesting to consider species with overlapping generations, where previously established individuals may block immigration into patches at carrying capacity (i.e., a priority effect; Atkins and Travis 2010). Priority effects could slow replacement of initial colonizers at the range edge, impeding genetic rescue and increasing the persistence of expansion load away from the edge. Some factors that could speed expansion beyond rates seen in our model and lead to increased expansion load include greater long-distance dispersal and increased fecundity. We found that simulations with fecundity halved (3.5 vs. 7; data not shown) exhibited much less expansion load and higher fitness in edge populations. Thus, increasing fecundity allows populations to persist at lower mean fitness and thus accumulate more load. Finally, other factors could result in complex, less predictable effects. Local adaptation likely involves multiple quantitative traits adapting over similar or different environmental gradients, with potentially complex interactions that merit future investigation.

A potentially key evolutionary component of range expansions not included in our model is the evolution of dispersal ability. Increased dispersal is always expected to evolve at expanding range margins (Hargreaves and Eckert 2014). Interestingly, increased dispersal is expected theoretically and found empirically even during expansion across environmental gradients or with expansion load alone (Henry et al. 2015). However, increased dispersal also steepens the perceived slope of a given environmental gradient, which can eventually slow or even temporarily halt range expansion until edge populations evolve to overcome initial maladaptation (e.g., Phillips 2012). Therefore, it is unclear how dispersal evolution would affect the results presented in our study.


Our results support those of previous studies finding that expansion load via the surfing of deleterious alleles reduces fitness in expanding populations (Peischl et al. 2013, 2015; Peischl and Excoffier 2015). We show this under biologically realistic conditions, bolstering evidence that allele surfing may indeed cause expansion load in nature. Our results are also in agreement with those of previous studies showing that on an environmental gradient, migration load reduces fitness in expanding populations (Kirkpatrick and Barton 1997; Bridle et al. 2010; Polechová and Barton 2015), although we did not see any cases of stable range limits as a result of local maladaptation or expansion load. We highlight the mechanism of interaction between local maladaptation and expansion load. Local maladaptation feeds back to reduce the speed of expansion and thus allows for less expansion load through reduced genetic drift and increased migration to the edge throughout the course of range expansion. Finally, we demonstrate that faster range expansion leads to a larger contribution of moderate- and large-effect deleterious alleles to expansion load. These contributions significantly advance theory on the genetics of range expansion toward meaningful predictions and interpretations for studies of natural populations.

We would like to thank D. Irwin, M. Kirkpatrick, S. Wang, and members of the Otto and Whitlock labs for feedback on early stages of this project. We also thank S. Peischl and an anonymous reviewer for insightful comments that improved the manuscript. Funding was provided by the Beaty Biodiversity Research Centre, Natural Sciences and Engineering Research Council (NSERC) Discovery Grant RGPIN-2016-03779 to M.C.W.; NSERC and Killam postdoctoral fellowships to N.P.S.; an NSERC Discovery Grant to A.L.A.; Swiss National Science Foundation (SNSF) Doc.mobility project P1SKP3_168393 to R.M.-D.; and SNSF grant PP00P3 144846 to F.G.; the Genome Canada Large Scale Applied Research Project Program supported G.L.C.

Appendix A Supplementary Results and Analyses

Table A1.

Fitness and expansion-speed results per scenario parameter combinations

   Fitness measures, core/edge
Dispersal kernel, environmental gradient (b)Genomic deleterious mutation rate (UD)Mean expansion speed ± 1 SETotal w¯Quantitative trait w¯zDeleterious loci w¯D
 .0.03.741 ± .005.932/.932.932/.9321/1
 .0375.02.314 ± .011.929/.399.929/.3991/1
 .375.0.756 ± .003.870/.199.870/.1991/1
 .0.13.559 ± .008.863/.748.933/.933.925/.802
 .0375.12.283 ± .011.859/.390.929/.461.924/.845
 .375.1.744 ± .002.804/.193.870/.217.924/.891
 .01.02.285 ± .014.496/.306.933/.945.532/.323
 .03751.01.991 ± .011.494/.268.929/.694.532/.387
 .3751.0.595 ± .002.463/.160.870/.366.532/.437
 .0.05.689 ± .010.933/.929.933/.9291/1
 .0375.02.695 ± .014.929/.369.929/.3691/1
 .375.0.736 ± .002.871/.183.871/.1831/1
 .0.15.355 ± .030.862/.770.933/.936.924/.823
 .0375.12.637 ± .019.859/.360.929/.420.924/.858
 .375.1.721 ± .002.803/.181.870/.206.923/.881
 .01.03.377 ± .045.498/.283.933/.938.534/.302
 .03751.02.285 ± .014.496/.249.929/.644.534/.388
 .3751.0.577 ± .004.463/.162.870/.367.532/.442
Figure A1.
Figure A1.

Range edge fitness predicts expansion speed, but more strongly in cases with strong environmental gradients. The speed of expansion increases more rapidly as fitness increases with decreasing steepness of the environmental gradient (black points and solid line; quadratic regression β1=−4.70, β2=9.15, P<.0001). The speed of expansion increases, but less drastically, with higher fitness from decreasing UD and reducing expansion load (open circles and dashed line; linear regression β=1.95, P<.0001).

Figure A2.
Figure A2.

Average core and edge fitness components for the quantitative trait and deleterious alleles in the case of the leptokurtic dispersal kernel (long-distance dispersal). Error bars indicate 95% confidence intervals, and points are offset for better visualization.

Figure A3.
Figure A3.

Recovery from expansion load, measured as wD over time from initial colonization to completion of simulations. Thin lines show all 20 replicate simulations. Dashed lines are for simulations with a kurtotic dispersal kernel (i.e., increased long-distance dispersal), while solid lines are for those with a Gaussian dispersal kernel. The mean for each set of 20 replicates is indicated by the thick line. Fitness was measured through time at the cross section in the center of the landscape (x=1,000). Each line begins at the time point when individuals were first recorded in the central cross section; thus, the starting point depends on expansion speed. The end point in all cases is after the landscape filled, which took longest on the steepest environmental gradient (green).

Figure A4.
Figure A4.

Proportion of expansion load due to loci fixed only at the range edge. Load for these loci is measured both at the range edge and in a segment of the range core with approximately the same total population size. Because we modeled only 1,000 loci subject to unconditionally deleterious loci, fixation for small-effect loci is inflated. Hence, in this analysis, we control for loci that are fixed both in the range core and at the range edge. We remove those loci from this analysis and calculate the load due to loci fixed only at the edge. Fixation is measured for the same population of individuals defined to be at the edge as in our fitness calculations. The amount of load due to fixed loci was much higher at the range edge because the alleles were at higher frequency in the edge populations. The proportion of load due to these fixed loci increased with decreasing steepness of environmental gradient and also increased with decreasing genome-wide deleterious mutation rate, UD. Both of these results make sense in light of our study’s findings that faster expansion contributes to increased expansion load. This shows that faster expansion also allows for increased fixation of deleterious alleles at the range edge.

Appendix B Supplementary Movies

Time Lapse of Average Fitness Components across the Landscape

Video B1.
Video B1.

Fitness across the landscape through time for each component of fitness, wD (deleterious mutations) and wz (quantitative trait), in the case of UD=0.0. Time is measured in generations relative to the end of the burn-in period of the simulation; that is, t=0 is the start of expansion. Each individual replicate (n=20) is shown, with the average across replicates indicated by the thick, darker points. Because different replicates have slightly different expansion speeds, the average at the range edge is noisy through space (video B1).

Video B2.
Video B2.

Fitness across the landscape through time for each component of fitness, wD (deleterious mutations) and wz (quantitative trait), in the case of UD=0.1. Time is measured in generations relative to the end of the burn-in period of the simulation; that is, t=0 is the start of expansion. Each individual replicate (n=20) is shown with the average across replicates indicated by the thick, darker points. Because different replicates have slightly different expansion speeds, the average at the range edge is noisy through space (video B2).

Video B3.
Video B3.

Fitness across the landscape through time for each component of fitness, wD (deleterious mutations) and wz (quantitative trait), in the case of UD=1.0. Time is measured in generations relative to the end of the burn-in period of the simulation; that is, t=0 is the start of expansion. Each individual replicate (n=20) is shown with the average across replicates indicated by the thick, darker points. Because different replicates have slightly different expansion speeds, the average at the range edge is noisy through space (video B3).

Time Lapse of Allele Frequencies across the Landscape

Video B4.
Video B4.

Frequency of deleterious alleles across the landscape throughout the course of a simulation. This simulation was selected as one where a large-effect allele fixes on the landscape (leptokurtic dispersal kernel, UD=0.1, b=0). Loci are colored by bins of homozygote effect size, as indicated in the key that appears in the video. Darker blues represent larger-effect loci. A single locus is colored in orange to represent one of the large-effect loci (s=0.0898) that fixed on part of the landscape during the course of expansion. Time is measured in generations relative to the end of the burn-in period of the simulation; that is, t=0 is the start of expansion. As time proceeds, it can be seen that more loci are fixed at the range edge as a result of allele surfing. Over time during expansion, recovery from fixation at the expanding front occurs because of immigration from behind the expanding front, as well as (minimally) from back-mutations, reducing the frequency of larger-effect loci (video B4).

Literature Cited