Consequences of Range Contractions and Range Shifts on Molecular Diversity

Abstract

Due to past and current climatic changes, range contractions and range shifts are essential stages in the history of a species. However, unlike range expansions, the molecular consequences of these processes have been little investigated. In order to fill this gap, we simulated patterns of molecular diversity within and between populations for various types of range contractions and range shifts. We show that range contractions tend to decrease genetic diversity as compared with population with stable ranges but quite counterintuitively fast range contractions preserve higher levels of diversity and induce lower levels of genetic differentiation among refuge areas than slow contractions. Contrastingly, fast range shifts lead to lower levels of diversity than slow range shifts. At odds with our expectations, we find that species actively migrating toward refuge areas can only preserve higher levels of diversity in refugia if the contraction is rapid. Under slow range contraction or slow range shift, active migration toward refugia lead to a larger loss of diversity as compared with scenarios with isotropic migration and may thus not be a good evolutionary strategy. These results suggest that the levels of diversity preserved after a climate change both within and between refuge areas will not only depend on the dispersal abilities of a species but also on the speed of the change. It also implies that a given episode of climatic change will impact differently species with different generation times.

molecular diversity, range contraction, range shift, range expansion, refugia, migration, climate change

Introduction

Climate changes constitute a crucial threat for the survival of many species (e.g., Easterling et al. 2000; Thomas et al. 2004). Range contractions or range shifts may occur as a consequence of temporal climatic fluctuations, depending on the geographical structure of the landscape, the duration of the climate change, or a species dispersal abilities (e.g., Hewitt and Yohe 1996; Walther et al. 2002; Parmesan and Yohe 2003), but one common response of animal species to climate change is migration toward more suitable regions.

Overall, these series of population expansions and contractions could play fundamental roles in the history of a species and in the emergence of new ones. During a contraction, a species’ range can be fragmented due to climatic heterogeneities and irregular geographic environments, leading to different subpopulations eventually residing in small isolated regions. Such regions are usually described as refuge areas (or refugia) (Stewart et al. 2010) and processes of population divergence follow, sometimes resulting in allopatric speciation events (Hewitt and Yohe 1996; Lister 2004; Jakob et al. 2007). On the other hand, range shift processes are also very frequent during climate changes (e.g., Truong et al. 2007; Tsang et al. 2008; Hoban et al. 2010; Tsai and Manos 2010). Consequently, many studies have tried to analyze and predict the migration behavior and distribution of range-shifting species in correlation with climate variables such as temperature, humidity, or precipitation (e.g., Peterson et al. 2002; Pearson and Dawson 2003; Schwartz et al. 2006). However, although most models were able to predict potential distributions based on current climates, different predictions were often obtained with different approaches under scenarios of climate change (see further details in Araujo et al. 2005; Elith et al. 2010). Some general predictions on species distributions empirically related to climate variables could nevertheless be extracted by using consensus models or by ensemble forecasting (Araujo et al. 2005; Araujo and New 2007).

Range shifts have received increased attention since the emergence of concerns about climatic changes. Using spatial and temporal models, Desai and Nelson (2005) calculated the critical velocity of a range shift above which populations become extinct. McInerny et al. (2009) simulated range shifts in a lattice under a colonization–extinction metapopulation model where climate change was modeled by a variation in extinction probabilities. They found higher levels of genetic diversity in the middle of the range under stable climates, while a larger diversity was detected in the receding border under slow range shifts. The authors explain this observation by a mismatch in the place of origin and the place of survival of extant lineages, which is the same in expectation under a stable climate. Later, Atkins and Travis (2010) analyzed the effects of local adaptation under climate change by a spatially explicit individual-based model. They found that higher dispersal levels or shorter climate change duration lead to lower rates of population extinction. But again, none of these studies looked at the effect of climate change and species dispersal abilities on genetic diversity, which can eventually determine the fate of the populations (Hewitt and Yohe 1996; Miller et al. 2006). Very recently, Cobben et al. (forthcoming) showed that genetic diversity was eroded during range shifts in three scenarios of increasing temperature under a metapopulation model. They found a more pronounced decline of diversity when the variation of temperature was faster. However, the role of dispersal abilities of species was not investigated in that study. Despite their extreme importance, it is curious that the genetic consequences of range contractions and range shifts have been generally much less investigated than range expansions (see Ray et al. 2003; Excoffier et al. 2009; Petit 2010), even though they should occur at least as frequently. To our knowledge, a single coalescent-based study addressed this question (Leblois et al. 2006), where the authors showed that an instantaneous habitat reduction can lead to genetic diversity loss, as measured by the number of alleles and heterozygosity. In addition, the authors showed that this loss can dramatically increase with the time spent in the refugia after the contraction. However, since only the case of an instantaneous habitat contraction was considered, a more realistic study seems necessary to check the generality of these conclusions.

Current patterns of genetic diversity are best understood by considering the genealogical processes arising under a given demographic scenario. Coalescent theory is a retrospective approach describing the genealogy of genes currently sampled in some population until their most recent common ancestor (MRCA). For neutral genes, this process only depends on the past demography of the population (Kingman 1982), and the genetic diversity of a sample can then be easily and quickly generated by adding mutations along the genealogy from the MRCA. Moreover, it is an extremely good approximation to classical forward population genetics processes (Wakeley 2008). It also provides intuitive explanations to complex phenomena. For example, after a demographic expansion, gene trees show long terminal branches (where most of mutations occur) because coalescent events are less likely in the currently large population than in the small ancestral populations (e.g., Hudson 1990). Conversely, a demographic contraction should lead to an accumulation of recent coalescent events due to small population sizes at the end of the contraction, resulting in gene trees with short terminal and long internal branches, and therefore to an overall decline in genetic diversity (Nei et al. 1975; Sousa et al. 2009; Peter et al. 2010). However, patterns of genetic diversity may change when considering spatial constraints and migration patterns, such as those occurring during range expansions (Wegmann et al. 2006; Barton 2008). For example, it is known that a range expansion looks like a large demographic expansion only if a large number of migrants are exchanged between neighboring demes (Ray et al. 2003; Excoffier 2004). Therefore, for a fixed migration rate, populations with large carrying capacities can show typical patterns associated to a demographic expansion, but populations with low carrying capacities having gone through a recent range expansion can look like stationary populations or even show signs of a demographic decline (for a review of the effect of range expansions on genetic diversity, see Excoffier et al. 2009). It seems therefore important to check if spatial and temporal range contractions would have the same effect as pure demographic contractions (i.e., bottlenecks Nei et al. 1975) or if migration patterns between demes and toward refuge areas would also affect genetic diversity, like in the case of spatial expansions.

In this paper, we have extended the SPLATCHE framework, initially developed to study the genetic consequences of habitat heterogeneity and range expansions in a spatial context (Currat et al. 2004; Ray et al. 2010), to deal with habitat contraction and anisotropic migration. We used SPLATCHE to simulate DNA sequence diversity under different scenarios of range contractions and range shifts. We report patterns of molecular diversity at different scales (within and between refuge areas) obtained as a function of the speed of these processes and migration patterns between neighboring demes.

Materials and Methods

All range expansion, range contraction, and range shift scenarios have been simulated with the program SPLATCHE 2 (Ray et al. 2010), which has been already used in several other studies (e.g., Francois et al. 2008; Schneider et al. 2010). Under this framework, genetic diversity is simulated in two distinct steps. The first step consists in a forward simulation of the demography (range expansion and range contraction) of a subdivided population where demes are arranged on a 2D stepping-stone lattice model (Kimura and Weiss 1964) and where either isotropic or anisotropic migration between neighboring demes can occur. The second step is a classical backward coalescent simulation (Kingman 1982), where the neutral genetic diversity of gene samples drawn from the population is generated. This molecular diversity only depends on the demographic information that has been recorded in the forward simulation step. We now describe the peculiarities of the forward simulations, considering several models of range contractions and range shifts. We then present how genetic diversity was estimated from simulated data.

Demographic Simulations

Following Ray et al. (2003), who first studied the effect of range expansion on genetic diversity, we performed simulations on a 2D lattice of 51 × 51 demes, where we included two southern isolated refuge areas of 5 × 5 demes that were separated by an empty area (see fig. 1). We arbitrarily chose to start the colonization of this world from the western refugium (at position <2;2>), and this location could be considered as a southern introduction place, the location of a speciation event, or the remains of a previous range contraction. Prior to the expansion phase, and in keeping with Ray et al. (2003), the size of the ancestral population was set to 100 individuals, thus assuming very limited initial genetic diversity. Each generation and in each deme, the life cycle begins by a logistic growth phase followed by an emigration phase where individuals from occupied demes migrate to neighboring demes at rate m (but this rate can vary over space and time, see below). The density of each deme is logistically regulated with an intrinsic growth rate r (set constant as r = 0.8, which is commonly found in nature [e.g., Sibly and Hone 2002]) and a carrying capacity K (i.e., maximum local density) that can be set to different values (K = 100 and K = 0 in favorable and unfavorable environments, respectively) over space and time during the range contraction or the range shift, in order to study its effect on final levels of diversity. We have also studied the sensitivity of our results to departures from assumed demographic parameters by considering alternative scenarios, with different sizes for the refugia (by considering refuge sizes of 2 × 2 and 5 × 25 demes), different growth rates (r = 0.6 and r = 1.0), and different carrying capacities in favorable environments (K = 50 and K = 200).

FIG. 1.

Simulated landscape and spatial processes. (A) The simulation landscape consists in a 2D lattice made up of 51 × 51 demes. The lattice is virtually divided into 23 horizontal layers of 2 × 51 demes, plus a 5-deme thick layer containing two refuge areas of size 5 × 5 demes. The four gray dots in each refuge represent the sampled demes. In panes B and C, black areas represent demes occupied at different times of the three different spatial processes, and white areas are empty. The timelines runs from left to right on both lines of panes B and C. The range contractions and range shifts shown in panes B and C, respectively, follow a range expansion after a homogenization time (see text and fig. 2). This range expansion always starts in the middle of the left refuge area. Note that range contractions and range shifts can proceed at different speeds as described in figure 2.

The expected number of emigrants sent by the i-th deme during the migration phase is noted N_it m, where N_it is the density of deme i at time t, and the actual number sent by each deme is drawn from a Poisson distribution with mean N_it m. Note that demes on the edge of the simulated lattice send an average of (N_it m)/4 migrants to each neighboring demes, implying that these edges act as partially absorbing boundaries. In contrast to previous models of migration where individuals are just replaced by new immigrants (e.g., Slatkin 1977), emigrants and immigrant counts can differ. The number of emigrants indeed depends on the local deme size, whereas the number of immigrants depends on the size of neighboring demes. Potential migration unbalance is taken care of by the population logistic regulation step.

The SPLATCHE program was modified to allow for anisotropic migration, where different emigration probabilities could be assigned to the four neighboring demes. In isotropic migration scenarios, we used emigration rates m₁ equal to 0.05 or 0.01 in each direction (south, north, east, west) corresponding to a total emigration rate of 0.2 and 0.04, respectively. Note, therefore, that in that case individuals are emigrating with equal probabilities to each neighboring deme. In scenarios with anisotropic migration, we used a large emigration rate m₂ = 0.2 or 0.1 toward a particular direction (north or south, see below) and a smaller emigration rate m₃ = 0.01 toward the other directions, which results in total emigration rates of 0.23 and 0.13, respectively. This anisotropic migration scheme was applied in several scenarios of range expansion, range contraction, and range shift.

Range Expansions

Range expansions were simulated with two objectives: 1) the colonization of a landscape before the simulation of a range contraction or a range shift and 2) setting a baseline scenario against which to compare the specific effects of range contractions and range shifts. In scenarios with range contractions, the range expansion spans the entire landscape, whereas only four layers are colonized in an expansion before a range shift (see fig. 1C and supplementary fig. S1, Supplementary Material online).

In order to prevent any sectoring effect happening during range expansions (Hallatschek et al. 2007; Excoffier and Ray 2008), we allowed for a homogenization period before implementing range contraction or shifts: these latter events occurred at least 5,000 generations after the end of the range expansion (which was relatively rapid as it took around 550 generations for scenarios with the smallest emigration rates). During range expansions, carrying capacities were kept constant (K = 100) over time. For proper comparisons with scenarios of range contraction and range shift under anisotropic migration, scenarios of pure range expansions were also simulated with similar anisotropic migration (see below) occurring on a layer in the northern edge without further contraction.

Range Contractions

As mentioned above, range contractions are simulated after a homogenization period following a range expansion that leads to the colonization of the whole world. Although other procedures could be envisioned, we have chosen to simulate a range contraction by making the available world progressively smaller and smaller with time. At fixed times, we make a horizontal layer (of size 2 × 51 demes and located on the northern edge of the world) uninhabitable by setting the carrying capacity of its demes to zero. A series of 23 consecutive range contraction will thus progressively make the colonized world shrink, until only two refuge areas remain (see fig. 1B). Different range contraction speeds were modeled by varying times between two successive episodes of contraction: i.e., 10, 50, 100, or 500 generations (fig. 2). We also studied the effect of allowing for different times spent in the refuge areas, called hereafter “refugial isolation time,” and samples were examined 10, 100, or 1,000 generations after the end of the contraction.

FIG. 2.

Schematic representation of the main types of range contractions that were simulated. In each case, after an initial phase of range expansion (green), the species range remains constant (blue) for some time allowing for the homogenization of genetic diversity among populations. Then, a phase of range contraction occurs (red), which is followed by a final phase of evolution in refuge areas (black). Genetic diversity is sampled at the end of the four phases. We also distinguish scenarios where the total simulation time varies according to the duration of the range contraction (A) and scenarios where the total simulation time is constant despite variable contraction times that are compensated by a proportional change in the homogenization time (B). Note that other scenarios allowing for variable evolution times in refuge areas are not shown here but would be represented by variable lengths of the black bar.

As the total simulated time is the sum of the expansion, homogenization, contraction, and refugial isolation times (see fig. 2), different total times were obtained when range contractions were simulated with different speeds (fig. 2A). We also studied scenarios with the same total simulated time but differing in expansion or contraction speed by varying the homogenization time accordingly (fig. 2B).

We evaluated the impact of the way species move toward refuge areas by implementing either isotropic or anisotropic migrations. In scenarios with isotropic migration, we assumed that individuals have equal probability to move to favorable (K = 100) or unfavorable (K = 0) environments. In scenarios with anisotropic migrations, we assumed that individuals living on the edge of the range would have the ability to “sense” their environment and would have a greater probability to move in the direction of the refuge areas, as if there was a gradient of environmental quality (even though we did not implement such a gradient). Individuals not living on the edge were assumed to be unable to sense the environmental gradient and migration remained isotropic in the range core. Note that anisotropic migration is only implemented in the layer located at the receding front of the range. Therefore, this modeling of directional migration can be considered as a way to simulate individual response to a forthcoming decrease in habitat quality by preferential migration to high-quality habitats.

Range Shifts

Like in range contractions, range shifts are simulated after a phase of range expansion and homogenization and involve a change in carrying capacities above and below the current habitable range, the size of which was arbitrarily set to four habitat layers (of total size of 8 × 51 demes). During a range shift, we distinguish the expanding from the receding front. The expanding front is a layer where carrying capacity has just been increased (generally to K = 100) and which has thus been newly colonized, whereas the receding front is a layer where deme carrying capacities will be set to zero in the next shift. We thus simulated a cycle of range shifts by moving the species range from the bottom to the top of the simulated world and back to the bottom in a series of 40 consecutive shifts (20 to the north followed by 20 to the south, see fig. 1C and supplementary fig. S1, Supplementary Material online). The speed of the range shift was varied in the same way as for range contractions. Such back and forth range movements are often observed in nature as seasonal movements (e.g., Dingle 1982), but they could be considered as the tracking of favorable environments during a climate warming followed by a colder period (deMenocal et al. 2000; Parmesan and Yohe 2003). In all cases, samples were taken 100 generations after the return of the range to southern refugia, to allow for the successful colonization of these refuge areas. Like in the case of range contractions, we studied cases where the total simulated time could vary or was kept constant by adjusting the homogenization time.

Both isotropic and anisotropic migration scenarios were also simulated during range shifts. However, in the anisotropic case, a larger emigration rate was only implemented in the layer at the receding front and in the direction of the expanding front (see supplementary fig. S1, Supplementary Material online; case B) to mimic a scenario where individuals try to escape from harsh conditions and to keep with the rest of the population.

Note that in our simulation framework, the speed of the range shift does not depend on the dispersal abilities of the species, but we consider that the speed is imposed by the environment (e.g., climate change), which implies that we model a “suitable range” shift where individuals can live and disperse.

Coalescent Simulations

Coalescent simulations were performed after the forward demographic simulations, by using the recorded information on migration rates and deme densities. Under each scenario, we reconstructed the genealogy of eight population samples, all of them consisting in 25 DNA sequences of 1,000 bp. Four demes were sampled per refuge area (as shown in fig. 1), which allows us to study molecular diversity at three different scales: 1) within demes, 2) within refuge areas, and 3) between refuge areas. Mutations were added on the genealogy starting from a randomly chosen sequence at the MRCA, with rate μ = 3.3 × 10⁻⁷ per generation and per base pair under a Jukes Cantor substitution model of evolution (Jukes and Cantor 1969), and without recombination. Note that this mutation rate was chosen to generate appreciable levels of molecular diversity and to prevent saturation or total lack of mutations in all scenarios. A total of 1,000 coalescent simulations were carried out for each demographic scenario.

Statistics Summarizing Genetic Diversity

Several summary statistics were computed to evaluate the genetic diversity of the simulated data sets using the program ARLEQUIN ver. 3.5 (Excoffier and Lischer 2010). In particular, we computed different indices of molecular diversity such as the number of alleles (k), heterozygosity (H), the number of segregating sites (S), the number of pairwise differences (π), and Tajima’s D (D) (Tajima 1989a), and this both at the deme and the refuge area level. We also measured genetic differentiation among groups with the F_ST and F_CT statistics computed under an AMOVA framework (Excoffier et al. 1992).

Results

The results are divided into two sections bearing, respectively, on the analysis of range contractions and on the analysis of range shifts. In each section, we describe the joint genetic consequences of different speed of climate change and different migration patterns of a species.

Effect of Range Contractions on Molecular Diversity

Our results show that the speed of the range contraction has important and somewhat counterintuitive effects on patterns of molecular diversity. Briefly (but see below), we find that 1) slow range contractions lead to a more severe reduction of genetic diversity than fast contractions and 2) active migration toward refuge area only leads to a better preservation of genetic diversity than isotropic migration for fast contractions but leads to even larger diversity loss if contractions are slow.

Range Contraction with Isotropic Migration

Globally, and as expected under a pure demographic contraction, range contractions decrease genetic diversity as compared with population with stable ranges. However, quite unexpectedly, we find that under isotropic migration, slow contractions result in much more reduced diversity than fast contractions (supplementary table S1, Supplementary Material online) for populations having high emigration rates (m₁ = 0.05). However, when emigration rates are low (m₁ = 0.01), these trends are not observed if the total simulated time varies according to the speed of the range contractions (see fig. 2A). This is because in that case more mutations occurred in the genealogies of the samples when the contraction is slower, as the total simulation time is also longer. Note that with a constant total simulation time, the decline in diversity is observed for any emigration rate (see supplementary table S4, Supplementary Material online and text below). Thus, when m₁ = 0.05, the average number of alleles $\stackrel{-}{k}$ strongly declines with the duration of the contraction, passing from $\stackrel{-}{k}$ = 29.7 ± 0.25 for a total contraction time of 230 generations to $\stackrel{-}{k}$ = 23.6 ± 0.22 for a total contraction time of 11,500 generations (supplementary table S1, Supplementary Material online). Note that both the number of segregating sites and the average number of pairwise differences increased with the duration of the contraction also because the total simulated evolutionary time increased in this case (see fig. 2A). Tajima's D values increase with the duration of the range contraction, passing from negative values after a rapid contraction, which is usually indicative of a demographic expansion, to positive values indicative of demographic contraction (Tajima 1989b) in scenarios with the slowest range contractions (supplementary table S1, Supplementary Material online). Therefore, signals of demographic decline evidenced by Tajima’s D are only visible here for the slowest contractions, and signals of old expansions are preserved by fast contractions. Slow range contractions also lead to higher levels of genetic differentiation between the two refuge areas. For example, with an isotropic emigration rate of 0.05 and a very slow range contraction of 11,500 generations, F_ST between refuges reaches 0.22, whereas it is only 0.13 for a fast contraction of 230 generations. However, as expected, low levels of migrations lead to very high levels of differentiation between refugia, irrespective of the speed of the contraction. Finally, note that levels of population differentiation within refuge areas do not depend much on the speed of the contraction (supplementary table S1, Supplementary Material online; F_ST levels when sampling only in the SW refugium).

Active Migrations Toward Refugia during a Range Contraction

In case of anisotropic migration, when individuals on the contraction front have a higher probability to migrate in the direction of the refuge areas, levels of genetic diversity are usually lower than in the case of isotropic migration when the contraction is slow (supplementary table S2, Supplementary Material online). In that case, the number of segregating sites is becoming even smaller with an increasing duration of the contraction, especially with high migration rates, despite the total simulation time being longer (fig. 2A). Tajima’s D is also more positive for slower contractions and thus gives stronger signals of demographic decline. F_ST estimates are found higher than with isotropic migration, indicative of stronger isolation between subpopulations. Interestingly, we find that the strategy of sending more migrants toward refuge areas only leads to high levels of genetic diversity in those refugia when contractions are rapid (supplementary table S2, Supplementary Material online). Otherwise, this strategy will always lead to a larger loss of diversity than if individuals migrate in random directions.

As could be expected, smaller anisotropic emigration rates toward the south during the range contraction produced results closer to those obtained under scenarios with isotropic migration (supplementary table S3, Supplementary Material online).

Effect of Refugial Isolation Time

As noticed previously by Leblois et al. (2006), the refugial isolation time deeply influences patterns of molecular diversity at both population and group levels: longer times in the refugia generally decrease levels of molecular diversity (fig. 3 and supplementary figs. S2 and Supplementary Data, Supplementary Material online) and the genetic differentiation between the refugia strongly increases (fig. 4 and supplementary figs. S2 and Supplementary Data, Supplementary Material online). Note that the effect of the speed of the contraction on diversity is erased by long refugial isolation times (e.g., >1,000 generations).

FIG. 3.

Average number of alleles (k) observed in refuge areas after range contractions of different speeds. Upper graphs: Total number of alleles observed among all pooled samples. Lower graphs: Number of alleles observed among the samples of the SW refugium (results for SE refugium are similar but not shown). (e) indicates cases where there has been no range contraction after the initial range expansion. Bar colors represent results obtained for different refugial isolation times: black, 10 generations; gray, 100 generations; white, 1,000 generations. The left graphs correspond to scenarios simulated with isotropic migration at rate m₁ = 0.05, whereas right graphs correspond to scenarios simulated under anisotropic migration (m₁ = 0.05, m₂ = 0.20, m₃ = 0.01). Error bars indicate 1.96 × standard error intervals obtained from our 1,000 simulations.

FIG. 4.

Average F_ST observed in refuge areas after a range contraction of different speeds. Upper graphs: F_ST values observed between SW and SE refuges. Lower graphs: F_ST values observed between samples within SW refugium (results for the SE refugium are similar but not shown). Note the zoomed y-axis on the lower graphs. (e) indicates cases where there has been no range contraction after the initial range expansion. Bar colors represent results obtained for different refugial isolation times: black, 10 generations; gray, 100 generations; white, 1,000 generations. The left graphs correspond to scenarios simulated with isotropic migration at rate m₁ = 0.05, whereas right graphs correspond to scenarios simulated under anisotropic migration (m₁ = 0.05, m₂ = 0.20, m₃ = 0.01). Error bars indicate 1.96 × standard error intervals obtained from our 1,000 simulations.

Controlling for Total Simulated Evolutionary Time

Because previous results could have been affected by differences in the total simulated time between slow or rapid contractions, we repeated our simulations by adjusting the homogenization time to have a similar total simulated time in all scenarios (fig. 2B). In that case, we find that slower contraction lead to more diversity loss for all summary statistics (indeed including the number of segregating sites and the average number of pairwise differences) and to larger extent of differentiation between refuge areas (supplementary tables S4 and Supplementary Data, Supplementary Material online).

Effect of Range Shift on Molecular Diversity

In range shift scenarios, shift speed and migration rates within the occupied range also have a large impact on resulting patterns of diversity. In our simulations (with different total times, like in fig. 2A), the fastest range shifts overall lead to the lowest levels of genetic diversity, in contrast to range contractions where lowest levels of diversity are observed for the slowest contractions (compare, e.g., figs. 3 and 5).

FIG. 5.

Average number of alleles observed in refuge areas after range shifts occurring at different speeds. Upper graphs: Total number of alleles observed among all pooled samples. Lower graphs: Number of different alleles observed among the samples of the SW refugium (results for the SE refugium are similar but not shown). (e) indicates cases where there has been no range shift after the initial range expansion. Bar colors represent results obtained for different migration rates: gray, m₁ = 0.05; white, m₁ = 0.01. The left graphs correspond to scenarios simulated with isotropic migration, whereas right graphs correspond to scenarios simulated under anisotropic migration (m₂ = 0.20, m₃ = 0.01). Note that no data are available when the range shift lasts only 400 generations in case of isotropic migration and a small migration rate m₁ = 0.01 due to population extinction (see text). Error bars indicate 1.96 × standard error intervals obtained from our 1,000 simulations.

Therefore, slow range shifts better preserve original genetic diversity but also lead to less differentiated populations within and between refugia (figs. 5 and 6, supplementary figs. S4 and Supplementary Data, Supplementary Material online). As expected, scenarios with low levels of gene flow between neighboring demes (m₁ = 0.01) lead to lower final diversity and higher population differentiation than scenarios with larger emigration rates (m₁ = 0.05) (figs. 5 and 6 and supplementary figs. S4 and Supplementary Data, Supplementary Material online). Migration is not only important for preserving genetic diversity, but it can also prevent population extinction, which was found to occur in scenarios of fast range shift with low isotropic emigration rates (supplementary fig. S6, Supplementary Material online).

FIG. 6.

Average F_ST observed in refuge areas after range shifts occurring at different speeds. Upper graphs: F_ST values observed between SW and SE refuges. Lower graphs: F_ST values observed between samples within SW refugium, SE refugium produced similar results (not showed) to the SW refugium. (e) indicates cases where there has been no range shift after the initial range expansion. Bar colors represent results obtained for different migration rates: gray, m₁ = 0.05; white, m₁ = 0.01. The left graphs correspond to scenarios simulated with isotropic migration at rate, whereas right graphs correspond to scenarios simulated under anisotropic migration (m₂ = 0.20, m₃ = 0.01). Note that no data are available when the range shift lasts only 400 generations in case of isotropic migration and a small migration rate m₁ = 0.01 due to population extinction (see text). Error bars indicate 1.96 × standard error intervals obtained from our 1,000 simulations.

The comparison of range shift scenarios with isotropic and anisotropic migration showed a somewhat surprising result. Unlike scenarios with isotropic migration where genetic diversity increases monotonously with slower range shifts (fig. 5, left column; supplementary fig. S4, Supplementary Material online), very slow shifts can lead to a decrease in diversity as compared with faster shifts in case of anisotropic migration (compare results for T = 20,000 and T = 4,000 generations in fig. 5, right column and in supplementary fig. S5, Supplementary Material online). We shall explain this phenomenon in the Discussion.

Finally, like in the case of range contractions, the effects reported above for range shifts remain qualitatively similar when the total simulated time is kept constant over scenarios (supplementary tables S6 and Supplementary Data, Supplementary Material online).

Effect of Alternative Demographic Parameters on Molecular Diversity

We have studied the influence of various growth rates (r = 0.6, r = 0.8, and r = 1.0), carrying capacities (K = 50, K = 100, and K = 200) and refuge area sizes (sizes 2 × 2, 5 × 5, and 5 × 25 demes) on the genetic diversity of both range contraction and range shift scenarios under isotropic or anisotropic migration. We find that different growth rates do not significantly affect the estimated pattern of genetic diversity within or between sampled demes (see supplementary figs. S8 and Supplementary Data, Supplementary Material online). However, we find that larger carrying capacities and larger refugia result in larger levels genetic diversity within demes and lower levels of differentiation between demes and refuge areas (see supplementary figs. S10–Supplementary Data, Supplementary Material online). However, despite these differences in overall diversity levels, the response of diversity to varying speeds of range contractions and range shifts remained qualitatively similar to what was described above for r = 0.8, K = 100 and refugia size of 5 × 5 demes.

Discussion

High Impact of Slow Contractions on Genetic Diversity

Our simulations show that the speed of contraction plays an important role on final levels of genetic diversity after a range contraction. Contrary to the naive view that strong and fast range contractions would have the highest impact on a species genetic diversity, our results rather show that fast range contractions better preserve initial levels of diversity and lead to little genetic differentiation between refuge areas. Contrastingly, slow range contractions lead to more reduced levels of diversity and stronger genetic differentiation between refugia. These results are in agreement with those of Leblois et al. (2006), who showed that an instantaneous contraction had little effects on the genetic diversity of remaining populations.

These results can be explained by considering the spatial and temporal distribution of coalescent events and by realizing that backward in time a range contraction corresponds to a range expansion. Therefore, with slow contraction, the area where gene lineages can potentially migrate backward in time will only grow very slowly and in that case it is likely that two genes might migrate to the same deme where they can coalesce. This event is much less likely in the case of a fast contraction because going backward in time, the range where genes can potentially migrate will increase very rapidly, and genes will have less opportunity to migrate to the same deme, and therefore, coalescence times will be longer. In other words, after a slow contraction, sampled individuals will be more closely related to each other than after a rapid contraction because they will have common ancestors living in the refuge area or very close to it. Thus, slow contractions should lead on average to trees with shorter terminal branches and smaller overall length than fast contractions and hence to less mutations and lower overall genetic diversity. With the same line of reasoning, levels of differentiation between refugia are explained by the fact that compared with two genes sampled from different refugia, two genes sampled in the same refuge area will have a relatively more recent common ancestor in case of slow contractions than when contraction is fast, which will lead to smaller F-statistics (Slatkin and Voelm 1991). The same effect explains why the refugial isolation time strongly reduces genetic diversity and increases the genetic differentiation among refugia: going backward in time, different genes sampled in the same refugium are likely to migrate to the same deme and hence coalesce within the refugia after the contraction. Note that a high genetic differentiation between refugia may promote speciation events as observed in real cases (e.g., Hewitt and Yohe 1996; Jakob et al. 2007).

Overall, we see that like pure population bottlenecks, range contractions reduce genetic diversity and that longer contractions like longer bottlenecks will have more effect (Nei et al. 1975). However, as was found in the case of range expansions, patterns of migrations between neighboring demes will modulate this diversity loss (see below).

Fast Range Shifts Have a Negative Impact on Genetic Diversity

In contrast to range contractions, our simulations show that fast range shifts lead to lower levels of diversity than slow range shifts. This is in keeping with results obtained by Cobben et al. (forthcoming) where faster temperature variations led to larger loss of diversity in shifting metapopulations. As can be seen from the distribution of coalescent events on supplementary figure S7 (Supplementary Material online), ancestral lineages tend to coalesce in narrow vertical bands located above the refuge areas for fast shifts, which implies that the effective range is more restricted in case of fast than slow shifts. Moreover, between two consecutive fast shifts, the ancestral lineages will have less time to leave the expanding front and will thus tend to remain on this expanding front instead of exploring the whole available range, hence inducing a faster rate of coalescent events, shorter trees, and less final diversity. In other words, diversity is also reduced in range shifts by recurrent founder effects on the expanding front (which actually does not exist in a range contraction) and populations have little time to recover from these founder effects if shifts are too fast. Note that in range shifts, the whole population is also at risk of extinction, if its dispersal abilities do not allow it to catch up with the fast moving receding front. In that case, individuals who do not follow fast range shifts are detached from the viable range (where K > 0) and die (see, e.g., Atkins and Travis 2010), which is indeed what we observed in case of very fast range shifts and low migration rates (see supplementary fig. S6, Supplementary Material online). Such population extinctions are in keeping with other studies on climatic change modeling (McInerny et al. 2009; Atkins and Travis 2010), and the observation that less resilient species are at higher risk of extinction (e.g., Burns et al. 2003; Berke et al. 2010).

Importance of Migration to Maintain Diversity

Although it is known that migration patterns may play an important role for the survival of species during climate change scenarios (Midgley et al. 2006), we are not aware of any other work that analyzed the effect of migration on genetic diversity after range contractions or range shifts. Our simulations showed that the migration abilities of a species are fundamental for maintaining its genetic diversity under both range contraction and range shift processes. As expected, species with higher migration rates always maintained higher levels of initial diversity under any range contraction or range shift scenario.

Active Migrations Toward Refuge Areas May be a Counter Productive Strategy

However, the beneficial aspect of migrations on genetic diversity has to be buffered by our finding that directional migration toward refuge areas (or higher quality habitat) is not necessarily a good strategy. Contrary to our prior expectation, which was that active migration toward refuge area would allow species to bring genetic diversity from the whole range into refuge areas, we see that in case of slow contractions, active migration toward the refugia leads to a more pronounced loss of diversity compared with scenarios with isotropic migration. This fact can be explained by an increased rate of coalescent events occurring on the receding front during the contraction. Indeed, the receding front, which is a source of emigrants in the forward process, becomes a sink for ancestral lineages when considering the process backward in time and genes will rapidly coalesce there, decreasing diversity. Therefore, contrary to our expectations in case of anisotropic migrations, genetic diversity present far away from the refuge areas is actually less efficiently brought to the refugia when contractions are slow because more diversity is lost during this contraction phase than when individuals move at random in all directions. Interestingly, lower levels of diversity are also observed in case of slow range shifts with anisotropic migration. In that case, ancestral lineages that could escape the founder effect occurring on the expanding front can freely roam the available range but are becoming attracted by the receding front if they wander just below it. This fatal attraction explains the drop in diversity observed for very slow range shifts (fig. 5). It is only in case of recent fast range contractions that active migrations toward refuge areas is beneficial and better preserves diversity than when individuals migrate in random directions on the receding edge (fig. 3). The effects of anisotropic migration during range shifts and contractions have been little studied so far, and further work would be necessary to understand their full effects. For instance, it could be envisioned that active migrations toward refuge areas would have a wider beneficial effect if they were distributed over the whole species range rather than restricted to the receding front or that some appropriate levels of directional migration could allow a species to preserve diversity whatever the speed of the contraction.

Robustness of Our Results to Varying Demographic Parameters

Our relatively limited sensitivity analysis showed that population growth has basically no effect on resulting patterns of diversity. However, diversity within samples increases with the size of the refuge area and with local deme densities, without modifying the dependence of diversity to the speed of environmental changes. Therefore, we believe that our results about the effects of range contractions and range shifts on diversity are likely to be qualitatively valid for a wider range of demographic parameter values.

Envisioning More Complex Scenarios

However, more complex demographic and ecological scenarios could have some influence on genetic diversity during range shifts and range contractions. For instance, local population dynamics where population growth varies across space and time or some form of nonreversible dispersal (density-dependent emigration and/or immigration) processes could alter the ability of mutations to surf to high frequencies in newly colonized habitats (Munkemuller et al. 2011) and, consequently, have an effect on final levels of diversity. Moreover long distance dispersal events, which are known to affect patterns of diversity after a range expansion (Nichols and Hewitt 1994; Ibrahim et al. 1996; Ray and Excoffier 2010), would certainly play a role during range contraction, most probably by buffering the negative effects of fast range contractions and shifts on diversity. Finally, we have only considered here a single round of expansion and contraction, whereas it is likely that many species have gone though several glacial cycles (e.g., Kropf et al. 2003; DeChaine and Martini 2004). We would therefore think that the phenomena we have described here would be amplified, but not qualitatively changed, by considering several consecutive cycles of expansions and contractions.

Evolutionary Implications

The fact that final levels of genetic diversity vary depending on the speed of contractions or range shifts implies that species with different generation times will be affected differently by the same episode of climatic change. Similarly, species with different dispersal abilities or migration strategies will loose various amounts of genetic diversity (Austerlitz et al. 2000). Our results thus suggest that generalist species with short generation times will suffer more from a given climatic change than those with longer generations, especially if they have limited dispersal abilities (Hamrick et al. 1992; Austerlitz et al. 2000). Conversely, specialist species that must track precisely their environment, and which should thus rather undergo range shifts than simple contractions, should be more affected by a given climatic change episode if they have long generation time and limited dispersal abilities. Moreover, a strategy consisting in migrating toward high-quality habitats should be more beneficial for species with long than with short generation time, as the best strategy for the latter seems to be migrating randomly. If one assumes that range contractions have occurred repeatedly in the evolution of most species and that life histories preserving genetic diversity are favored (e.g., Booy et al. 2000), one could postulate that long-lived organisms should have preferentially developed migrating strategies allowing them to track the most suitable environment.

We want to thank Matthieu Foll and Gerald Heckel for constructive comments made during this study and other CMPG laboratory members for helpful discussions. We also thank two anonymous reviewers for very insightful comments and suggestions. This work was supported by the EU/FP7 grant N° 226852 (SCALES).

Author notes

Associate editor: Ryan Hernandez