Modelling the effect of habitat fragmentation on range expansion in a butterfly

Robert J Wilson, Zoe G Davies, Chris D Thomas

Abstract

There is an increasing need for conservation programmes to make quantitative predictions of biodiversity responses to changed environments. Such predictions will be particularly important to promote species recovery in fragmented landscapes, and to understand and facilitate distribution responses to climate change. Here, we model expansion rates of a test species (a rare butterfly, Hesperia comma) in five landscapes over 18 years (generations), using a metapopulation model (the incidence function model). Expansion rates increased with the area, quality and proximity of habitat patches available for colonization, with predicted expansion rates closely matching observed rates in test landscapes. Habitat fragmentation constrained expansion, but in a predictable way, suggesting that it will prove feasible both to understand variation in expansion rates and to develop conservation programmes to increase rates of range expansion in such species.

Keywords:

1. Introduction

Biodiversity is not static, but responds dynamically both to natural and anthropogenic changes to the environment, with population and distribution responses often lagging behind environmental changes that have already taken place (Tilman et al. 1994; Brooks & Balmford 1996; Helm et al. 2006). Small remnant populations and metapopulations (networks of local populations connected by dispersal; Hanski 1999) are prone to extinction through chance events, so there is no guarantee that protecting remnant populations will conserve threatened species successfully, even in the absence of further habitat degradation (Hanski & Ovaskainen 2002; Bulman et al. 2007). Therefore, conservation programmes must attempt to bring about genuine species recovery. In addition, species responding to climate change may be threatened if they are unable to expand their ranges through heavily modified landscapes to reach new areas of suitable climate space (Thomas et al. 2004), and there is evidence that many species are failing to keep up (Warren et al. 2001; Menéndez et al. 2006). To manage the distributions of rare and threatened species efficiently, given limited time and resources, quantitative means are needed to predict how distributions of suitable habitats will affect rates of expansion at species range margins.

Metapopulation models represent one means of modelling distribution change for species using fragmented or patchy habitat networks (Hanski 1999). In a metapopulation, the extinction rate of local populations is higher in small or low-quality habitat patches, and the colonization rate of empty patches is low for patches that are isolated from existing populations. Therefore, metapopulation models predict that the probability of metapopulation persistence (Hanski & Ovaskainen 2000), the proportion of patches occupied at equilibrium (Hanski et al. 1995) and the rate of range expansion through fragmented habitat (Thomas et al. 2001a) will be greater for landscapes with a greater overall area and connectivity of habitat. Expansion beyond current species range margins may be inhibited by landscapes in which a lack of habitat leads to low rates of colonization, and to high rates of extinction in populations that become temporarily established (Bascompte 2003; Opdam & Wascher 2004). Conservation programmes need to convert these general expectations into quantitative predictions. Metapopulation models make many simplifications, but if their predictions are nevertheless accurate they have the advantage of being relatively simple to parametrize, and therefore practical in a way that models requiring detailed life-history data may not be, given that such data are unavailable for virtually all rare species of conservation concern.

In this paper, we apply a metapopulation model to the recent dynamics of a rare butterfly that has expanded its range, as habitat availability and local population sizes have increased in association with habitat management and climate change (Davies et al. 2005, 2006). We estimate parameters for the model in a training landscape where the species had a relatively stable distribution between 1982 and 2002, and use the parameters to model species dynamics in four independent testing landscapes where it expanded its distribution over the same period. Comparison of observed and modelled changes shows that metapopulation simulations can lead to accurate predictions of the relative likelihood of occupancy of different habitat patches, and of rates of expansion in landscapes differing in habitat quality and configuration.

2. Material and methods

(a) Study system

The silver-spotted skipper butterfly Hesperia comma L. is confined in England to species-rich calcareous grasslands, where it lays eggs exclusively on short tufts (less than 10 cm) of sheep's fescue grass Festuca ovina (Thomas et al. 1986). Hesperia comma experienced severe habitat loss in the mid-twentieth century as a result of agricultural intensification (which has destroyed over 80% of Britain's calcareous grassland since 1940; Asher et al. 2001), and vegetation succession following the abandonment of livestock grazing and decline of rabbit populations through myxomatosis. A full survey of H. comma populations in 1982 showed that the species had declined to fewer than 70 populations covering a total area of only 2.1 km2, in eight habitat networks in southern England (Thomas et al. 1986). Here, we focus on the five networks that contained all but three of the 1982 H. comma populations: the North Downs in Surrey and Kent, the South Downs in Sussex and Hampshire, and the Chiltern Hills (figure 1).

Figure 1

The distribution of H. comma and its habitat in five networks in Southeast England. Symbols show 2×2 km grid squares containing suitable patches: black squares, H. comma populations present 1982 and 2000; grey squares, patches colonized 1982–2000; pluses, population introductions 1982–2000; open circles, all patches vacant 1982–2000; open triangle, two populations extinct 1982–2000. More than one patch may be present in each grid square. Solid line shows English coast.

Since 1982, habitat restoration has taken place for H. comma as a result of recovering rabbit populations, and active conservation management including the re-establishment of livestock grazing in many areas where H. comma was absent in 1982. Habitat availability and local population sizes may also have increased as a consequence of climate change, apparently because a wider range of microhabitats now satisfy the species' thermal requirements than in 1982 (Davies et al. 2006). Despite these improving conditions, H. comma has only spread within habitat networks where it was already present, and more than 90 per cent of colonizations have been within 10 km of 1982 populations (maximum 29 km; Davies et al. 2005). Here, we test whether differences in range expansion rates in different regions are governed principally by variation in habitat availability.

(b) Habitat and species distribution

Distribution surveys were carried out in all five networks in 2000, and in Surrey and Sussex in 2002. All suitable chalk grassland containing tufts of F. ovina shorter than 10 cm was surveyed within 30 km of the 1982 H. comma distribution, and within 10–15 km of newly colonized sites found in 2000 (see Davies et al. 2005). Habitat patches were defined as areas of suitable grassland bounded by continuous woodland or scrub barriers, or by at least 25 m of unsuitable grassland (Thomas & Jones 1993; Hill et al. 1996). Hundred metre grid coordinates, area, altitude, aspect, slope, shelter and vegetation composition (including percentage cover of F. ovina growing in turf less than 10 cm tall) were recorded for each habitat patch. The H. comma flight period in each network was determined using weekly transect counts at a reference population. All habitat patches were searched for adults and eggs during the regional flight period, with at least one visit towards the end of the flight period to search for easily found eggs in patches where H. comma was not found on the first visit.

The maximum distance from 1982 populations to any colonized patch in 18 years was just below 30 km, so we estimate network habitat availability as the area of suitable habitat in 2000 within 30 km of 1982 populations. Habitat quality for each patch was defined as the percentage cover of F. ovina growing in vegetation less than 10 cm tall. To test effects of variation in habitat quality on range expansion, patch area was adjusted by dividing the percentage cover of suitable F. ovina in each patch by the average percentage cover of suitable F. ovina in all patches and then multiplying this term by patch area.

(c) Metapopulation modelling

Metapopulation models provide an appropriate framework for H. comma's regional dynamics because (i) H. comma breeds in clearly defined, localized habitat patches (Thomas et al. 1986), (ii) most populations are small and at some risk of extinction, with occasional extinctions in 1982–2000 despite generally increasing population sizes (Thomas & Jones 1993; Davies et al. 2005), (iii) low rates of dispersal link H. comma populations (Hill et al. 1996) and colonization rates decrease with distance from established populations (Thomas & Jones 1993; Davies et al. 2005) and (iv) the probability of habitat patch occupancy increases with patch area and connectivity to other populations of the species (Thomas et al. 1992).

We used the incidence function model (IFM; Hanski 1994, 1999), a stochastic patch occupancy model based on the assumptions that extinction rate declines with habitat patch area, and that colonization rate increases with patch connectivity. Connectivity (Si) for a habitat patch i is defined as Embedded Image (Moilanen & Nieminen 2002). Parameter α is the slope of a negative exponential dispersal kernel, where the proportion of per generation dispersal over distance d km or greater corresponds to expαd; dij is the distance to patch i from each occupied source patch j (where ij); and Aj is the area (ha) of each patch j. Source patch emigration rate scales with patch area to the power b, which was set to 0.5 to account for the tendency of per capita emigration to be greater from smaller habitat patches in this and other species (Hill et al. 1996; Moilanen & Nieminen 2002).

In the IFM, annual colonization probability of patch i is Embedded Image, and annual extinction probability, Embedded Image. Patches with higher connectivity (Si) are more likely to be colonized, and patches with larger area (Ai) are less likely to go extinct; 1−Ci simulates the rescue effect by the instantaneous re-colonization of patches that would otherwise go extinct. Extinction risk is unity where minimum patch area A0=e1/x. Thus, A0, e and x scale the relationship between patch area and extinction, and y determines the relationship between connectivity and colonization probability. Based on the assumptions that extinction is area dependent and colonization is connectivity dependent, IFM parameters e, y, x and α can be estimated from snapshots of patch occupancy using maximum-likelihood estimation (Moilanen 1999, 2000). Software for parameter estimation and metapopulation simulation was obtained from http://www.helsinki.fi/science/metapop/.

For this study, we estimated IFM parameters e, y, x and α using the distributions of occupied and vacant H. comma habitat patches in Surrey in 2000 and 2002. The distribution of the species has remained far more stable in Surrey than in other networks (Thomas & Jones 1993; Davies et al. 2005), but nevertheless overall habitat availability has increased (Davies et al. 2006), so we estimated metapopulation parameters using the two distribution snapshots closest together in time, when habitat area is unlikely to have changed markedly (2000: 74 occupied patches, 32 vacant; 2002: 78 occupied, 28 vacant). We used standard parameter estimation techniques, assuming that the two snapshots of the Surrey distribution were representative of a stochastic steady state (Moilanen 1999). We set the minimum area at which populations could survive from one year to the next as A0=0.02 ha, based on field observations of occupied habitats, and assumed zero regional stochasticity. Parameters estimated using the Monte Carlo Markov Chain procedure, with 1000 function evaluations in initiation and 4000 function evaluations in estimation, were x=0.415, y=8.885, e=0.198 and α=0.405.

Parameter estimates from Surrey were used to run 100 simulations of 18 years each, starting with H. comma's distribution in each network in 1982, and assuming that all habitat mapped as suitable in 2000 was available for colonization from 1982 onwards. This simplification is acceptable because much of the habitat was already available but unoccupied in the 1980s and early 1990s (Thomas & Jones 1993), although habitat management and climate change would have increased habitat availability over the 18-year study period (Thomas et al. 2001a; Davies et al. 2005, 2006). We assumed no colonization from outside the networks, which were separated by more than 30 km from other H. comma populations. For each network, the likelihood that each patch would be colonized was estimated from 100 simulations, using (i) total patch area and (ii) habitat-adjusted area, based on F. ovina density. The fit between observed and modelled colonizations for each network was analysed using the area under the curve (AUC) for receiver operating curves in SPSS v. 15.0 for Windows. Range expansion rate in each network (for observed and modelled data) is summarized as the mean distance to each patch colonized by 2000 from the nearest patch that was occupied in 1982. For different patches, this colonization distance may be measured from different 1982 populations, although in practice always from populations in the same network. This approach should lead to relatively conservative but realistic measures of range expansion, since patches near one extreme of a network are likely to be colonized from nearby populations rather than from elsewhere in the network.

(d) Analytical estimation of expansion rate

We also estimated the potential rate of H. comma's expansion in the absence of habitat barriers using an analytical method (Van den Bosch et al. 1990; Lensink 1997). The annual velocity of population expansion (Cexp) can be predicted based on the average distance from hatching that an individual produces offspring (σ), the mean age of first reproduction (μ), the expected number of offspring produced per individual in its lifetime (R0), the standard deviation of the first age of reproduction (v) and the kurtosis of dispersal distances (γ), using the equation (Lensink 1997)Embedded ImageFor H. comma, an insect with annual non-overlapping generations, μ=1 and v=0. Assuming a negative exponential distribution of dispersal distances, σ=1/α (the slope of the negative exponential distribution) and γ=0 (kurtosis for a normal distribution). Therefore, using the same negative exponential distribution with slope α=0.405 as in IFM modelling, σ=2.47 km. R0 has been estimated for H. comma as 1.5 (Hanski & Thomas 1994). Hence velocity of population expansion was estimated as Cexp=2.47*Embedded Image=1.42 km per year.

3. Results

A total of 732 habitat patches were identified in the five habitat networks in 2000. Out of 732 habitat patches, 243 contained H. comma populations: 62 populations surviving from 1982, 174 apparent colonizations and 7 resulting from reintroductions (patches colonized as a result of reintroductions occurred beyond the range of natural expansion and were excluded from analysis). Two 1982 populations (both in the Chilterns) had suffered extinctions and 487 patches remained unoccupied (table 1). The two extinct populations were predicted to suffer extinction in more than 50 per cent of IFM simulations, using either observed or habitat-adjusted areas, unlike any other populations occupied in 1982.

View this table:
Table 1

Occupancy, colonization and habitat availability in five Hesperia comma habitat networks, and the amount of variation explained in patch colonization patterns by metapopulation modelling. (Notes: significance levels: **p<0.01; ***p<0.001.)

Mean distance to patches colonized by 2000 from the nearest population in 1982 ranged from 1.3 km in Surrey (maximum 6.7 km) to 5.8 km in Sussex (maximum 28.8 km; table 1). The distance expanded in each network was positively but not significantly related to the total area of suitable habitat in 2000 (Pearson r2=0.54, n=5 networks, p=0.16). However, the average per cent cover of the host plant F. ovina growing in suitable conditions was significantly less in the Kent network than in each of the other four networks (Mann–Whitney U tests, p<0.001 for comparisons between Kent and patches in each other network). When habitat quality was taken into account using habitat-adjusted areas (table 1), the distance expanded in each network was significantly related to habitat availability (r2=0.81, n=5, p=0.04).

Observed colonizations in each network after 18 butterfly generations (years) were significantly related to IFM predictions of each patch's relative colonization probability, with network AUC values of 0.80–0.91. Adjusting areas for habitat quality improved the AUC for patch colonization in Kent, but otherwise did not improve predictions of patch colonization patterns (table 1). The proportion of habitat patches colonized declined with increasing distance from 1982 populations, but the declines with distance were not gradual, and differed markedly among networks (figure 2). The model's quantitative predictions of the proportion of patches colonized in 2 km distance intervals from 1982 populations were generally accurate, again with the exception of Kent (figure 2). IFM predicted relative expansion rates in each network well (figure 3; r2=0.95, n=5, p=0.001; regression through the origin, B=0.84±s.e. 0.09), but significantly overestimated expansion in the low-quality Kent network (median modelled average distance expanded =5.7 km, 95% CIs 5.2–6.6 km; observed=3.6 km). Simulations using habitat-adjusted areas no longer overestimated expansion distance in Kent (median=4.0 km, 95% CIs 3.2–4.8 km; figures 2f and 3). Habitat-adjusted simulations were significantly related to observed network expansion rates (r2=0.94, n=5, p=0.001; B=0.93±s.e. 0.12), but slightly underestimated the expansion in the Chilterns (median=1.6 km, 95% CIs 0.8–3.5 km; observed=3.8 km).

Figure 2

Proportion of habitat patches colonized by H. comma over 18 years in 2 km distance intervals from the nearest populations in 1982. Full network simulations: (a) Surrey, (b) Sussex, (c) Chilterns, (d) Kent and (e) Hampshire. Habitat-adjusted simulation: (f) Kent. Filled triangles show observed proportion of patches occupied in 2000. Solid lines show median of 100 simulations; dashed lines show 2.5th and 97.5th percentiles. Open triangles in (d) and (f) show patches occupied following introductions in 1997.

Figure 3

Average distance of patches colonized by H. comma in 2000 from patches occupied in 1982. Modelled expansions show median of 100 simulations, error bars from 2.5th to 97.5th percentiles. (a) Total area simulations and (b) habitat-adjusted simulations. Square, Surrey; triangle, Hampshire; diamond, Chilterns; large circle, Sussex; small circle, Kent. Solid diagonal lines show observed equal to modelled expansion.

Using the analytical procedure assuming that the only constraints to range expansion were H. comma's dispersal pattern and intrinsic rate of increase, we estimated that the species would expand radially at a rate of 1.4 km per year, leading to an expansion of 25.5 km after 18 years. Observed maximum patch colonization distances (from the nearest population in 1982) were less than 10 km in all networks except the Chilterns (max=16.6 km) and Sussex (max=28.8 km; table 1); and beyond 10 km only a small fraction of available habitat was colonized (figure 2).

4. Discussion

This study used a metapopulation model incorporating little of the detailed ecological and behavioural attributes of a test species beyond knowledge of the distribution of its habitat and populations in the landscape. Nonetheless, metapopulation parameters estimated from H. comma's occupancy patterns in one habitat network (Surrey) led to accurate predictions for four other networks both of distance expanded (figures 2 and 3) and the relative likelihoods that individual patches would be colonized (table 1). When combined with simple information on habitat quality (see Thomas et al. 2001b), the model captured sufficient essence of the species range expansion that it could estimate species recovery accurately and independently in different landscapes, and identify a landscape where low habitat quality constrained expansion rate.

Relatively simple metapopulation models can predict rates and patterns of range expansion well because they focus on the critical processes of colonization and extinction. If such models are well parametrized, then their predictions may be as accurate as those based on more complex approaches such as individual-based models (Ovaskainen & Hanski 2004). Predictions are also dependent on a sufficiently accurate measurement of the landscape-scale distribution of suitable habitat. Our assumption that all habitat identified as suitable in 2000 was available for colonization in 1982 was a simplification, and the incorporation of dynamic patterns of habitat availability in metapopulation models could be vital for modelling changes to species distributions (Keymer et al. 2000), especially in the context of climate change.

In 2000, H. comma was found to occupy 21 km2 of suitable habitat (10 times its area of occupancy 18 years earlier), in more than 250 discrete populations (an increase from fewer than 70 in 1982) (Davies et al. 2005). However, although 174 patches were colonized in the five habitat networks, and only two populations went extinct, in each network more patches remained unoccupied in 2000 than had been colonized since 1982 (table 1). An analytical estimation of range expansion by H. comma through continuous habitat greatly overestimated the average colonization rate of the species in 18 years, whereas metapopulation modelling based on stepping-stone colonization through fragmented habitat networks led to rather accurate predictions of range expansion rates. Thus, limited habitat availability constrained expansion rates substantially, over and above the constraints imposed by the intrinsic dispersal capacity and rate of increase of the species. The fact that habitat fragmentation constrained recovery, even in an expanding species with a landscape-scale conservation programme in place, implies that fragmented landscapes are likely to present an almost insurmountable barrier to the distributional responses of many species to climate change. Continuing the metapopulation simulations for H. comma for 100 years did not lead to complete patch occupancy in any network: many isolated habitat patches at the margins of each network were occupied in fewer than 50 per cent of simulations, and the section of the North Downs between Surrey and Kent (figure 1) was never predicted to be colonized naturally (R. J. Wilson 2008, unpublished data).

The results shed light on why, despite pronounced increases in population size on monitored sites (Davies et al. 2005), H. comma has failed to expand northwards in recent years to anything like the extent achieved by many British butterflies that have less specialized requirements (Warren et al. 2001). In this respect, H. comma shares similarities with many butterflies (Menéndez et al. 2006) and other taxa (Hickling et al. 2006) that are expanding their distributions at their poleward range boundaries, but are apparently spreading much more slowly than would be expected if they were keeping track with climate change. Previous modelling studies of the effect of habitat quantity on range expansion rates have concentrated on species with widespread habitats (e.g. Hill et al. 2001). Although these studies show that habitat availability influences expansion rate, such species are fairly atypical: habitats for most species occupy only a very small percentage of the landscape (Cowley et al. 1999). These localized species face a far greater challenge responding to climate change.

Modelling the responses of species distributions to habitat availability can identify landscapes where conservation is either most needed or most likely to improve species status (Huxel & Hastings 1999). Such an approach is potentially applicable to a wide range of systems, assessing the potential for natural expansion in response to habitat restoration, the likelihood that species introductions and reintroductions will succeed and spread, and the potential responses of species to climate warming (Opdam & Wascher 2004). Habitat fragmentation represents a major impediment to the expansion of species range boundaries as the climate warms (Warren et al. 2001; Travis 2003). Many species may only be able to change their distributions in response to climate change in landscapes where there is sufficient density of habitat to allow expansion (Peterson et al. 2002). Approaches similar to those described here can be used to model rates of expansion in different landscapes, to predict those species that will be able to shift their distributions in response to climate change, and those that will require active conservation management to do so.

Acknowledgments

We thank H. Burton, K. Ericson, P. Ewin, R. Fox, S. Glencross, A. Goodhand, S. Hanna, D. Hoare, C. Holloway, R. Leaper, J. Mellings and A. Moilanen for their assistance with fieldwork and analysis. Two anonymous referees provided helpful comments on an earlier draft. Funding was provided by the UK Natural Environment Research Council, EU TMR FRAGLAND and English Nature.

Footnotes

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

    • Received May 27, 2008.
    • Accepted September 29, 2008.

References

View Abstract