Dynamics of range margins for metapopulations under climate change

B.J Anderson, H.R Akçakaya, M.B Araújo, D.A Fordham, E Martinez-Meyer, W Thuiller, B.W Brook

Abstract

We link spatially explicit climate change predictions to a dynamic metapopulation model. Predictions of species' responses to climate change, incorporating metapopulation dynamics and elements of dispersal, allow us to explore the range margin dynamics for two lagomorphs of conservation concern. Although the lagomorphs have very different distribution patterns, shifts at the edge of the range were more pronounced than shifts in the overall metapopulation. For Romerolagus diazi (volcano rabbit), the lower elevation range limit shifted upslope by approximately 700 m. This reduced the area occupied by the metapopulation, as the mountain peak currently lacks suitable vegetation. For Lepus timidus (European mountain hare), we modelled the British metapopulation. Increasing the dispersive estimate caused the metapopulation to shift faster on the northern range margin (leading edge). By contrast, it caused the metapopulation to respond to climate change slower, rather than faster, on the southern range margin (trailing edge). The differential responses of the leading and trailing range margins and the relative sensitivity of range limits to climate change compared with that of the metapopulation centroid have important implications for where conservation monitoring should be targeted. Our study demonstrates the importance and possibility of moving from simple bioclimatic envelope models to second-generation models that incorporate both dynamic climate change and metapopulation dynamics.

Keywords:

1. Introduction

Shifts in the distributions of species with climate change have now been documented for many species (Warren et al. 2001; Hickling et al. 2005; Root et al. 2005; Parmesan 2006; Rosenzweig et al. 2008) and many more are expected to shift with future climate change (Thomas et al. 2004). To date, studies examining the likely effects of climate change on species distributions, and implicitly the range limits of these species, have used bioclimatic envelope models (Pearson & Dawson 2003; Ohlemüller et al. 2006; Thuiller et al. 2008). These studies generally use empirical models (e.g. Thuiller 2003; Segurado & Araújo 2004; Elith et al. 2006) to infer a statistical relationship between the current distribution of the species and the climate where it occurs. This relationship is used as a basis to project the future distribution of the species based on climate forcing predictions from individual runs or ensembles of global climate models (e.g. www.cgd.ucar.edu/cas/wigley/magicc/index.html).

However, one limitation of this procedure is that it does not account for population or metapopulation dynamics that determine species distributions, population structure and extinction risk at local scales (Thuiller et al. 2008), thus it may give misleading results when used to project the extinction rates (Akçakaya et al. 2006). Accounting for dynamics in demography and dispersal is particularly important at the limits of species ranges, where small population size, isolation, release from density-dependent intraspecific competition and limitations to inter-population movement are likely to play non-trivial roles in determining not just the presence but also the persistence of a species (Hanski & Thomas 1994; Keitt et al. 2001).

Keith et al. (2008) have recently used a coupled bioclimatic envelope/metapopulation approach to investigate the effects of various factors on the population viability of South African fynbos plants under stable and changing climate scenarios. We extend this approach to investigate the influence of metapopulation dynamics and dispersal on the range limits, as well as extinction risk, under stable and changing climate change scenarios. We investigate how populations at the expanding and trailing range margins differ from each other and from the core of the species distribution. Theory and limited empirical data (e.g. Caughley et al. 1988) suggest that shifts in population abundance and demographic rates along the ‘edge of the range’ should be some of the first and most sensitive signs of a broader species response to environmental change. As a case study, we used two lagomorphs (Romerolagus diazi and Lepus timidus) in different regions (Mexico and Great Britain) with different distribution patterns.

In response to climate change, species can shift, adapt or die (Foden et al. 2007; Brook et al. 2008). Shifts can be either polewards (e.g. for the Northern Hemisphere; Hill et al. 1999; Parmesan et al. 1999) or upslope (e.g. Shoo et al. 2005; Merrill et al. 2008). Here, we investigate the predicted poleward shift of L. timidus and the upslope shift of R. diazi, and the interactions between the population-demographic and climate response at their range margins.

2. Material and methods

We chose two lagomorph species with contrasting geographic patterns as study systems: the range-restricted R. diazi (volcano rabbit) and the locally abundant British metapopulation of the relatively widespread L. timidus (European mountain hare; see the electronic supplementary material for details). Spatially explicit metapopulation models were developed for R. diazi and L. timidus separately, using a three-step process.

(a) Environmental suitability models

An environmental suitability map was constructed for each species using a combination of GIS map layers based on habitat (land use and/or land cover, i.e. vegetation; see the electronic supplementary material for details) and climate (based on the species's modelled climatic envelope; see the electronic supplementary material for details of statistical fitting procedure). The values for these relative suitability layers—ranging between 0 and 1—were multiplied together to create a combined environmental suitability map. A different climatic map was produced for each year between 2010 and 2100, based on the A2 HadCM3 projections, using linear interpolation (see the electronic supplementary material for details), allowing us to include a dynamic annual change in climate suitability. The land use/land cover map was static.

(b) Metapopulation model

Spatially explicit pre-breeding stage-based matrix models were constructed for each species in RAMAS GIS (Akçakaya & Root 2005) and parametrized life-history data were sourced from previous demographic studies and abundance time series. The models were female-based and incorporated demographic and environmental stochasticity, dispersal and density dependence. The environmental suitability layer (see above) was used to determine the spatial structure of the metapopulation and, more specifically, the distribution of suitable patches in each year, occurring above a minimum threshold value based on the distribution of observed occurrences and suitability values (see the electronic supplementary material for details). The grid-cell clustering algorithm is described in Akçakaya & Root (2005) and used a neighbourhood distance of 3 km.

The metapopulation-demographic model is linked to the environmental suitability layer (climate suitability and land-cover suitability combined) through the carrying capacity (K). This was calculated for each patch summing across all grid cells within the patch. To calculate K, suitability (land cover×climate) was multiplied by an estimate of maximum density in optimal habitat, where land-cover suitability equalled 1 above a threshold value (50% cover for R. diazi and 70% for L. timidus) and 0 otherwise. Optimum density for L. timidus was set at 28 km−2, which was the value required to yield a total K of 150 000 females—a value estimated in a recent census of L. timidus in Great Britain (JNCC 2007). In a similar way to the mountain hare, the optimum density for the volcano rabbit was derived iteratively so that the model correctly predicted the population size of one of the largest populations (El Pelado), which was estimated as 6488 by Velázquez (1994). For both species, a habitat patch was considered too small to be suitable if it could not support a population size of at least 50 females—a rule of thumb value below which inbreeding depression and Allee effects (not modelled) often become important determinants of extinction risk (Akçakaya & Root 2005).

(c) Simulations

Population viability and species distributions were modelled under two scenarios—climate change and a reference baseline (stable climate). Initial population size and total initial carrying capacity were the same for each scenario. The climate change scenario was modelled by integrating a different climate envelope map for each year between 2010 and 2100. Consequently, species distribution, abundance and range limits were driven by demographic processes, climate change and the interaction between these. The stable climate scenario was modelled by keeping the climate envelope map for the first year (2010) static throughout the simulation. Thus, under this scenario, only demographic processes could cause changes in species distribution, abundance and range limits. For L. timidus, we modelled three dispersal rates (high, medium and low) for each scenario. This was done by altering the dispersal constant (b; see the electronic supplementary material).

All simulations were based on 10 000 replicates and run over a 90-year period (i.e. 2010–2100). The population viability was assessed using estimated female abundance, patch occupancy and cumulative probability of decline (Akçakaya & Root 2005).

(d) Range limits

The mean elevation and the elevational range limits of R. diazi were calculated as the average, minimum and maximum elevation of all the grid cells over the area occupied by all extant populations. Historic distribution was predicted based on records obtained from four scientific collections (see the electronic supplementary material for details). Core area was calculated as the total area of cells belonging to a patch for which all eight neighbours also belong to the same patch.

The northern and southern range limits for L. timidus were calculated based on a weighted mean of the latitudes of the most northern/southern 10 per cent of the metapopulation. Weights were the average population abundance of each patch in each year, and latitude was taken from the geographic centre of the patch. The total population centroid was calculated as the weighted mean latitude of the whole metapopulation.

3. Results

(a) Metapopulation dynamics

Romerolagus diazi and Lepus timidus show distinct differences when climate change is included in the suitability map (that which determined the spatial structure of the metapopulation model). Climate change amplifies the decline in female abundance and the cumulative probability of decline for R. diazi (figure 1), but only to a minor extent for L. timidus (figure 2). The number of occupied patches is also reduced somewhat with climate change for L. timidus (figure 2b), but for R. diazi the number of occupied patches increases from 2 to 4 by mid-century and then drops back to a mean of less than 2 (figure 1b). Closer inspection of the patch structure shows that this apparent increase in the number of patches is due to the existing larger patches fragmenting into smaller ones, rather than unoccupied patches becoming occupied (results in the electronic supplementary material).

Figure 1

Predicted status of the R. diazi metapopulation: (a) female abundance, (b) number of occupied patches and (c) cumulative probability of decline to fewer than 1000 females. For each diagnostic, two scenarios are presented: dynamic climate change (black line) and a static climate (grey line). Error bars, based on 10 000 stochastic model iterations, are presented for every 10th year only with the two scenarios offset by 1 year for clarity.

Figure 2

Predicted status of the L. timidus metapopulation within Great Britain: (a) female abundance, (b) number of occupied patches and (c) cumulative probability of declining to less than 10% of the initial abundance. For each diagnostic, two scenarios are presented: dynamic climate change (black line) and a static climate (grey line). Error bars, based on 10 000 stochastic model iterations, are presented for every 10th year only with the two scenarios offset by 1 year for clarity.

(b) Range limits

Metapopulation dynamic changes at the range margins are more pronounced, for both species, than shifts in the overall metapopulation.

When climate change is introduced to the suitability map, the lower elevational range limit of the R. diazi metapopulation increases steadily upslope, from approximately 2300 m to approximately 3000 m, which is just below the average elevation of 3150 m for the initial year (2010) and above the average elevation of 2800 m for the ‘historic’ distribution (figure 3). At the same time, both the total range area and the core area decrease steadily from 2010 onwards. The maximum elevation remains constant for both scenarios over time as R. diazi is already at its upper elevation range limit. Thus the overall effect under climate change was a range retraction. Under the stable climate scenario, the average, minimum and maximum altitudes did not change.

Figure 3

Predicted changes for the R. diazi metapopulation in: (a) total area (solid lines) and core area (dashed lines) of all occupied patches; and (b) mean (dashed lines) and minimum (solid lines) elevation. Two scenarios are shown: the dynamic climate change scenario (black) and the stable climate scenario (grey). Isolated points (his) represent an estimate of the historical area and elevation, respectively.

Under the climate change scenario, the L. timidus metapopulation shifts north over time (figure 4). However, this northerly shift is not equal throughout the distribution. The northern range margin shifts north almost thrice as far as the centroid of the metapopulation, and twice as far as the southern range margin. When combined with different dispersal estimates, the range margins respond quite differently to each other and to the centroid. Although the northern margin shifts north fastest under higher dispersal (figure 4a,d,g), the movement rate of the centroid is relatively unchanged (figure 4b,e,h) and the southern range margin shifts slower and therefore remains further south than the lower dispersal scenarios. The reason for this is apparent under the stable climate scenarios. With low dispersal, the southern range margin (figure 4c) remains relatively stable. But by contrast, with high dispersal, the southern range margin actually shifts south (figure 4i) and the northern range margin shifts north and then stabilizes. This is due to a dispersal-driven expansion, in both latitudinal directions, as movement from occupied to vacant patches allows a greater proportion of the total metapopulation to be colonized.

Figure 4

Predicted shifts in latitude for L. timidus under three different dispersal scenarios: (ac) low, (df) medium and (gi) high. Trends are shown for: (a, d, g) the northern range limit (weighted mean latitude of the northern 10% of the population); (b, e, h) the core centroid of the population (weighted mean latitude of whole population); (c, f, i) the southern range limit (weighted mean latitude of the southern 10% of the population). For each range margin and the population as a whole, two scenarios are presented: dynamic climate change (black line) and a static climate (grey line). Values represent deviations from initial values. Dot-dashed line represents 2010 values.

4. Discussion

Poleward and upslope movements of species with climate change have been documented in many different taxa and regions (e.g. Parmesan et al. 1999; Shoo et al. 2005). Many of these studies show advances or invasions (but see Wilson et al. 2005; Franco et al. 2006; Merrill et al. 2008). Climate envelope approaches are able to predict shifts in climate space in both leading and trailing edges (e.g. Thomas et al. 2004), but are unable to predict species persistence when the trailing edge shifts away from its present location (Akçakaya et al. 2006). Here, by directly linking spatially explicit predicted changes in climate to a metapopulation model, including demographic determinants of viability and dispersal, we are able to comment on both range limits, i.e. the leading and trailing edge of the distribution.

As a simplification, we used the same stage matrix (survivals and fecundities, and their variability) for the two species. When a model is density-dependent (as all of ours are), the details of the stage matrix may have some influence, but these details are not nearly as important as the intrinsic population growth rate, Rmax, and the variability in the vital rates. We used exactly the same Rmax for the two species, based on a weighted average of several lagomorph populations. The population dynamics that result from the two models are different, because the spatial structure, bioclimatic relationships, population sizes (or carrying capacities) and their distribution among patches, and the change in spatial structure (because of climate change) are all quite different for the two species.

Incorporating dynamic climate change into the R. diazi metapopulation had little effect in the first few decades, followed by a steep increase in the cumulative probability of decline. The resulting metapopulation is much smaller, and occupies a few isolated and very small patches at the end of the simulation (figure 1b) compared with the static reference. This means that the metapopulation will be much more susceptible to Allee effects, stochastic catastrophic events and changes in land use. None of these factors are considered in the current models but all of them are a very real risk to the continued persistence of R. diazi. However, Carroll (2007) has shown that although static climate change had the largest effect on lynx and marten declines, trapping and logging had an additive detrimental impact. The current distribution of R. diazi is a fraction of the historical distribution (especially in terms of ‘core areas’ of patches; see figure 3a), largely because of the changes in land use. Its range limits were considerably affected by the inclusion of dynamic climate change in the modelling (figure 3b). The top of the highest volcano the rabbit inhabits is 5452 m and the other volcano is 5222 m. Climatically, it would therefore have some room to move upwards. However, according to the static vegetation cover map we used, suitable habitat does not currently reach to the very top, so any movement upslope, as predicted by our models, corresponds to a loss of range. Although this may be an overly conservative assumption, vegetation may be constrained by other factors such as soil type and moisture availability, or shifts in suitable vegetation may simply lag behind climate change. Merrill et al. (2008) have shown a similar pattern of contraction at the lower elevation (attributed to climate factors), without a corresponding expansion at the upper elevation (attributed to biotic factors—lack of host plant) in the Sierra de Guadarrama (Spain) populations of Aporia crataegi.

In contrast to Keith et al. (2008), where the population viability of South African fynbos plants was relatively unaffected by dispersal, the response of the range limits of L. timidus to climate change was relatively sensitive to the degree of dispersal included in the model (figure 4). Although we included different estimates of dispersal as a sensitivity analysis, both empirical (e.g. Simmons & Thomas 2004) and theoretical studies have shown that dispersal ability may evolve differentially across the core-range margin gradient, increasing at the expanding range margin and decreasing where habitat is most fragmented (Hughes et al. 2007). With a high dispersal estimate, the northern range limit (figure 4g) expanded faster and further than with the low dispersal estimate (figure 4a,d), but the absolute expansion by the year 2100 is no different when compared with medium dispersal (figure 4d). However, L. timidus already occupies a fairly northern location within Great Britain. The shift rate of the metapopulation centroid was largely unaffected by the dispersal estimate used.

It is the trailing edge of species' distributions that is theoretically most affected by the population dynamics, changing more slowly for longer-lived and sessile species (Hampe & Petit 2005; Foden et al. 2007). Changes at this range margin are driven by a combination of movement and local extinction, whereas changes in the leading edge are largely dispersal-driven. Increasing the dispersive estimate of the species causes the metapopulation to respond slower, rather than faster, to the climate change on the southern range margin or trailing edge. This is probably because greater dispersal distances allow southern patches to be connected to the more northerly patches for a longer period.

There are many more instances in the literature of climate-induced shifts in the leading range margin of species, partly because it is more difficult to identify changes in the trailing edge of distributions with climate change owing to artefacts of sampling and recorder biases (Thomas et al. 2006; but see Wilson et al. 2005; Franco et al. 2006; Merrill et al. 2008). However, here, under a simulated situation of perfect information, we demonstrate that this lack of evidence might be due to genuine biological lags, even in shorter-lived, highly mobile species such as the mountain hare.

One of the pressures on the mountain hare appears to be competition from the European brown hare (Lepus europaeus), which is mediated by land cover, climate and interactions between these factors; the European brown hare is more competitive in open or agricultural land, especially in milder climates (Hewson 1991; Thulin 2003; Jansson & Pehrson 2007). Failure to include biotic interactions might further mask the effects of climate change on populations (Araújo & Luoto 2007) and their range limits, especially at the trailing edge, where competition from invading species is likely to have a detrimental effect. In addition, there is much uncertainty within our analysis (e.g. population parameters; climate projections; local adaptation; systematic differences in the dispersal ability; small pockets of climatically suitable habitat, which may exist at finer resolution). For the sake of simplicity and clarity, these are essentially ignored in the current study. However, when making management decisions, one needs to consider these sources of uncertainty, and future studies should attempt to improve our understanding of the range of possible outcomes and provide an estimate of the likelihood of these outcomes.

Coupling spatially explicit climate change with a dynamic metapopulation model is a substantial step forward in more realistically predicting the response of biodiversity to climate change. As we have shown, including this detail can reveal changes at the range margins of retreating and advancing geographic distributions that are overlooked by other methods. Further advances will come from simultaneously projecting changes in vegetation cover (animal habitat) and animal population dynamics in response to climate change, and, additionally, in simulating the multi-species effects of biotic interactions (e.g. competition or predation) on the persistence of metapopulations.

Acknowledgments

Thanks to David Roy (Centre for Ecology and Hydrology/Biological Records Centre) for help with the British distribution data regarding the European mountain hare and Yvonne Collingham (University of Durham) for help with the British climate data. We thank Andrew Smith for providing information on the volcano rabbit. Thanks to Ralf Ohlemüller, Aldina Franco, Chris Thomas, Georgina Mace, Steve Willis, Richard Pearson and Kevin Gaston for their useful discussion. B.J.A. was funded by UKPopNet (NERC R8-H12-01) and English Nature, and M.B.A. is funded by the Spanish Ministry of Science and Innovation (CGL2008-01198-E). M.B.A. and W.T. also received support from the EC FP6 MACIS (Minimisation of and Adaptation to Climate change Impacts on Biodiversity no. 044399) and EcoChange (Challenges in Assessing and Forecasting Biodiversity and Ecosystem Changes in Europe, no. 066866 GOCE) projects. This study is the result of two workshops on ‘Modelling species extinctions risk under climate change’ that took place at the Centre for Population Biology, Imperial College, London, and at the National Museum of Natural Sciences, CSIC, Madrid (the latter workshop was funded by the EC FP6 MACIS project). This work complies with UK law. B.J.A. obtained the L. timidus GB data and did the GIS, final figures, climate suitability modelling and lead writing. H.R.A., B.W.B. and D.A.F. did the RAMAS metapopulation modelling. B.W.B. conceived the manuscript. M.B.A. provided data on L. timidus. E.M.-M. provided data and information on R. diazi and generated the historic distribution model. W.T. provided the European climate data and assisted with the climate suitability modelling. All authors assisted with and commented on the manuscript.

Footnotes

  • One contribution of 17 to a Special Issue ‘Geographic range limits of species’.

    • Received September 29, 2008.
    • Accepted December 9, 2008.

References

View Abstract