Fragmentation and degradation of habitats are major contributors to the endangerment of species (Brooks et al. 2002, Reed 2004, Watling & Donnelly 2006). The impact of fragmentation depends mainly on two factors. Firstly, the physical configuration of the landscape structure, so called structural fragmentation (Tischendorf & Fahrig 2000). Secondly, how individuals perceive their surroundings and their ability to disperse, so called functional fragmentation. Fragmentation may lead to extinction, or to populations being prone to extinction, due to decreasing genetic variation and inbreeding depression, and/or stochastic ecological events (Gutzwiller 2002, Stacey & Taper 1992). If the deterministic threats leading to fragmentation can be stopped or reversed, a species may find opportunities to disperse between patches and reduce the risk of extinction (Landweber & Dobson 1999, Clobert et al. 2001). However, in order to stop or reverse deterministic threats, it is necessary to understand the dynamics of the population in question and to estimate which capacity is high enough to sustain a threatened species at the landscape level. A sustainable capacity is defined as a landscape's capacity to sustain a metapopulation (Hanski & Ovaskainen 2000).

A regional population consisting of several subpopulations (due to a patchy distribution of suitable habitats) connected by migration, and in which regional dynamics are driven by local extinction and recolonisation events, is characterised as a metapopulation (Hanski & Gilpin 1997). Although there are several factors influencing turnover in a metapopulation, e.g. birth and death rates, patch quality and predation, interconnection between subpopulations is a crucial factor for the existence of a metapopulation (Schumaker 1996). One of the more widely-used techniques in modelling metapopulation dynamics within real landscapes is the incidence function model (IFM; Hanski et al. 1996, Sjögren-Gulve & Hanski 2000, Hokit et al. 2001). In contrast to Levins' (1969) simple metapopulation model, IFM is a stochastic, spatially realistic, occupancy model describing the probability of occupancy of each single habitat patch (Hanski 1994, Hanski et al. 1996, Hanski & Gilpin 1997). It is possible to apply this model with few parameters; the information required is presence/absence 'snapshot' data from a set of habitat patches, which are defined according to patch area and pair-wise distances between patches. With knowledge of these parameters one can simulate and manipulate a patch network for a species (Hanski 1994).

The consequence of fragmentation on original habitat is an increase in that habitat's edge component and therefore an increase in the extent of edge effect (Primack 1998). The distance that edge effects extend into patches is highly variable (Laurence & Yensen 1991), and this may have consequences when unit area is used in estimations of population size and dynamics. Generally, the area of a patch is considered as linearly correlated with population size or species diversity (Connor & McCoy 1979, Martin & McComb 2003). In IFM, the extinction probability is assumed to be proportional to local population size which in turn is proportional to the local patch area (Hanski 1994). Although the basic area of patches is of importance, influence from the geometric shape of patches may be equally or more important in determining potential incidence of a species in a particular patch (Burrough 1981, Krummel et al. 1987, Helzer & Jelinski 1999, J. Sahlsten unpubl. data). The shape in terms of perimeter*area ratio may also be important when considering suitability of patches to sustain a viable population (Collinge & Palmer 2002), especially with management goals that consider connectivity and flow of genes, organisms and energy (Milne 1988, Jorge & Garcia 1997). If the shape of patches has an impact on the capacity of patches to sustain a local population size, it should also have an impact on the probability of extinction. Hence, basic area measurements might underestimate the impact of changes in a landscape, or perhaps overestimate the capacity of a certain patch network to sustain a metapopulation. Considering capacity in addition to occupancy levels is also important, because occupancy levels may be proportionally higher in a landscape although it has lower capacity, which may lead to overestimation of population viability in a landscape. One applicable method to use in dealing with regional dynamics in fragmented landscapes is metapopulation capacity (Hanski & Ovaskainen 2000). This method is commonly used in conservation planning or management strategies to rank different patch networks in terms of their capacity to sustain a viable metapopulation. However, only basic area measurements are currently used in these applications which can result in a biased estimation of population viability in a specific landscape and thereby risking misguided recommendations.

In this study we examine the effects of four different area-related scenarios, basic area and area reduced by varying amounts of edge zone, on occupancy level and capacity estimates using metapopulation theory. To alter the area*perimeter ratio and recalculate the effective area in our data we apply a GIS (Geographic Information System) method and reduce patch area by edge zones at different magnitudes by using negative buffering. By doing this, perimeter lengths, including hollows within patches, are automatically taken into account. Two patches with equal area but different perimeters are reduced proportionally to this difference, e.g. smooth vs convoluted perimeter. The study species is the hazel grouse Bonasa bonasia, a habitat specialist with poor dispersal ability and therefore likely to be strongly affected by habitat fragmentation (Andrén 1994, With & Crist 1995). Using GIS remote sensing, a technique ideal for monitoring large and remote areas, and the ecological niche factor analysis (ENFA) of Hirzel et al. (2002a,b) and Hirzel & Arlettaz (2003), a network of suitable patches for this species was identified (J. Sahlsten, unpubl. data). Occupancy level and capacity were estimated for this predefined patch network using our four different area-related scenarios. The results from these simulations indicate that the perimeter-area related size measure does not affect proportion of occupancy. On the other hand, the landscape seems to be sensitive to changes in terms of its capacity to support a viable metapopulation. We conclude from our results that perimeter-area related measures of patch size combined with capacity could be a more important measure for the estimation of population dynamics and impact of landscape changes compared to basic area measurement and occupancy levels. The implication of our results suggests that common and game species viability may be overestimated and consequently overexploited. In terms of predicted population declines for rare species, our results suggest that some species may be included in a lower endangerment category than they actually should be in the Red List of Threatened Species (IUCN), due to overestimations where basic area measurements have been used.

Material and methods

Incidence function model

In order to fit a patchily-distributed species to an IFM some assumptions are necessary (Hanski et al. 1996, Hanski & Gilpin 1997). All patches in the system are assumed to have a significant probability of going extinct. The dynamics of local patches are assumed to be much faster in comparison with the whole metapopulation and therefore local dynamics can be ignored. The patches in the network are assumed to be interconnected by migration. Furthermore, because this is an area where the hazel grouse has been established for thousands of years and there has been no major alteration of patch configuration in recent times, it can be assumed that populations of hazel grouse have reached a steady state of colonisation and extinction events. Given these assumptions, it was possible to model hazel grouse dynamics in the area as a metapopulation. The patch network and census results were used to set up the initial configuration of the patches and initial patch occupancies in the study area (Fig. 2).

The IFM is derived from a linear first-order Markov chain of presence or absence of a species in a habitat patch (Hanski 1994, Hanski et al. 1996). The unknown parameters within the model can be estimated with few parameters by applying presence/absence data from one survey (Hanski 1994, Hanski et al. 1996). Furthermore, it is also possible to make predictions of the impact on metapopulation dynamics from varying area-related measures of patch structure. The mathematical structure of the model involves the probability that a species, or even an individual, manages to colonise an empty patch and, at the other end of the scale, the probability of a patch going extinct. Thus the equation includes both probability of colonisation (C_i) and extinction (E_i) of the patch i. Using these terms, it is possible to calculate the incidence (J_i) of a species in patch i.

Incidence J_i of patch i is first defined in terms of extinction (E_i). The probability of extinction is related to the area (A_i) of the patch, and hence it can be used as a variable to determine extinction risk according to equations 2 and 3.

Thus, probability of extinction is indirectly determined by patch area (A_i) and by the environmental stochasticity (x). Secondly, incidence is defined by colonisation (C_i) rates. The colonisation rate is in turn affected by migration. Implying that rate of colonisation (C_i) can be derived from migration rate of the patch i (M_i).

where M_i is migration and y is an estimated constant that determines how fast colonisation probability reaches unity with increasing number of immigrants. Migration may also be expressed in terms of connectivity between patches (S_i). This connectivity is determined by factors such as presence of the species in a specific patch (p_i) and by the extinction probability (e). Further features that determine the connectivity are the distance between two patches (d_ij) and Alpha (α), which is a contant setting the survival rate for migrating individuals. Finally, the area of a patch (A_j) is important in terms of frequency of migration and thereby the area affects connectivity as well.

where

and p_j is one for presence and zero for absence in the patch from which a potential migrant comes. Alpha (α) is a constant setting the survival rate or inverse of dispersal distance ability for migrants over the Euclidian distance between patch i and patch j (d_ij). Area (A_j) reflects the emigration frequency from patch j to patch i. If equation 6 is substituted into equation 4, it gives a colonisation rate according to equation 7. ß is a constant that stands for a number of components which remain unknown in practice. This new equation contains a constant describing the ability of a species to colonise a patch (y′).

where y′ = y/ß describes the colonisation ability of the focal species. Lower case y′ indicates a good coloniser. Substituting expression 2 and 7 in equation 1 gives incidence (J) of the species in patch i:

To summarise, the data required to apply IFM are patch areas (A_i) and their spatial locations, in order to calculate the pair wise distances (d_ij), and the presence/absence of the species in the patches in the year of the survey (p_i). The remaining parameters y, e and x are estimated using equation 8. The fourth parameter, α, is species specific and fixed after estimation; the other parameters are estimated in relation to α. If the knowledge about dispersal patterns is insufficient, it is possible to estimate α from mark-recapture (Hanski 1994). In this study we set α to 0.25, which reflects a dispersal distance of 4 km. This dispersal distance was chosen for two reasons. Firstly, this is an average post-juvenile dispersal distance recently reported in the French Alps (Montadert & Leonárd 2006). Secondly, the patch configuration has an average nearest neighbour distance between patches of 1,629.2 m ± 597.7 S.D. (J. Sahlsten, unpubl. data). Thus, in order to sustain a metapopulation in this area, we assumed that hazel grouse have to disperse more than the 800 m reported by Swenson & Danielsen (1995).

There are several possible ways to create an IFM and, in our study, we modified a script developed earlier by Oksanen (1994) to fit the data we gathered, and applied it within the statistical software R (R Development Core Team 2005; Appendix I). By using appropriate tools in ArcMap (ESRI 2007), patch areas and nearest neighbour distance were obtained. Although the mean nearest neighbour distance is limited to landscapes of similar extent and resolution for comparative studies (Hargis et al. 1998), it can provide information of spacing between patches in a single landscape. Presence/absence data for hazel grouse were collected by a species census. However, it was necessary to estimate the other parameters (x, y and e) by approximation, with maximum likelihood regression incorporated in the R-script.

Reduced area

One basic assumption in the IFM is that the size of a subpopulation is correlated to the size of the patch, and extinction probability is correlated to the size of the population (Kindvall & Ahlén 1992, Hanski et al. 1996). The area alone could be a deceptive measure because some part of the area may be eliminated as suitable patch due to an edge effect. It has been shown that abiotic, biotic and ecological edge effect occurs at least 50 m from the edge towards the interior of a patch (Gates & Gysel 1978, Andrén & Angelstam 1988, Matlack 1993). In addition, it has been found that kill-sites of radio marked hazel grouse were closer to edges than to random sites (Swenson 1991b). Thus, in order to reduce patch area with the amount of potentially unsuitable area, three different edge zones of 50-, 75- and 100-m were applied to the original patch areas. As the proportion of edge-to-area increases, the quality of a patch in terms of available area is assumed to be reduced. Patches with a high perimeter value or high perforation level will consist of more edge/area than a patch with low perimeter or perforation values. Thus, if two patches have equal area but one has a higher perimeter value, the one with higher perimeter value will be reduced relatively more in terms of usable area compared to a patch with a lower perimeter. Consequently, the amount of area actually available will be determined by its perimeter to area relationship. Some of the original patches became < 20 ha after removing the edge area and were consequently removed from the patch network since they no longer were considered large enough to sustain a hazel grouse home range over the long term. Simulations of the IFM with 200 time increments were made with the basic area measurement and with the three reductions of area in terms of perimeter-area relationship. In order to control for variation of occupancy prediction by the model, simulations of two hundred time increments were looped two hundred times and the mean occupancy was extracted from the simulation result.

Metapopulation capacity is the sum of individual patch contributions which sustain a metapopulation in a landscape, and it can be used to predict consequences of degradation or loss of habitats on metapopulation dynamics (Hanski & Ovaskainen 2000). Close to equilibrium, a system can be described by a set of linear equations which may have special exponential solutions. Growth rates of these special solutions are called the eigenvalues, and the largest is called the leading eigenvalue. Capacity of a metapopulation is technically the leading eigenvalue (λ_M) of a landscape matrix with the elements:

The first component is reflecting the expected longevity of patch i, which depends on patch area (A_i) and some scaling of extinction (x) together with immigration rate (e). The second part A_j^em f(d_ij) corresponds to exp(-αd_ij)A_j, which is the connectivity component from equation 6, and it gives the rate at which patch i is colonised by patch j. However, since the absolute values of λ_M is dependent on units of measurement for area and distance, it is not possible to conclude if a particular value is high or low. On the other hand, with the assumptions of extinction rate = e/A_i and colonisation rate = c Σ_j≠i exp(-αd_ij)A_jp_j(t), different patch networks could be ranked according to their capacity. Hence, we ranked the capacity of the different patch networks to sustain a viable metapopulation via metapopulation capacity and manipulation of area measurement.

Results

Discussion

The number of occupied patches was lower in the patch network reduced by a 100-m edge zone compared to the basic patch network. This indicates that there is a risk that population sizes will be overestimated and consequent recommendations for management or conservation plans might be misguided. Further implications are that, although some patches are large enough to sustain viable local populations, it may be devastating for the regional population if patch area is overestimated. If one examines the proportion of patches occupied, there was no difference in occupancy between the area measurements. This result is an effect of concentrating suitable patches due to removal of initially empty patches that became too small to sustain a long-term hazel grouse home range after edge zone reduction. This implies that the amount of suitable habitat will also be overestimated if only the basic area of patches is considered.

This reasoning is supported by the fact that the patch network with reduced area shows a lower capacity to sustain a metapopulation. As a consequence, the habitats under consideration will be more sensitive to landscape changes than is indicated by basic area measurements. This higher sensitivity will have an important impact on planning strategies for forestry and conservation programmes. Additionally, our results indicated that many small changes in the landscape will have a cumulative impact that increases perimeter*area ratio or perforates patches until a threshold value is reached where the regional population collapses. Although we saw a clear difference in capacity between our different patch networks, which allowed us to rank them in terms of their ability to sustain a viable metapopulation, the differences in capacity could be even larger than the results indicate. This is because we used a fixed value for α in our simulations, when in theory α should increase as the amount of edge increases. However, further study of the relationship between dispersal and α are needed to determine how α should vary in the model.

Incidence function models have so far mostly been applied to model-species with rapid generation times, small patches and short inter-patch distances, such as butterflies, frogs and insects (Hanski et al. 1996, Biedermann 2000, Vos et al. 2000). In our study, we have attempted to apply the model to hazel grouse, a species which, like many other bird and mammal species, lacks the traits characteristic of model-species, which makes it more difficult to obtain the data needed for modeling. However, it is important to apply these models and methods to non-model species, even though it may be difficult to obtain data about colonisation-extinction events or when generation time spans are prohibitive, because these types of species are often of concern in conservation. It has, however, been shown that IFM predicted the population dynamics well for the American pika Ochotona princes using occupancy data (Moilanen et al. 1998), but there is still a need for further development of methods and models to estimate and predict population dynamics. In every study of a species distribution, it is essential to collect data about the species' preferred habitat structure. ENFA is a method that makes this possible with limited knowledge of the species ecology and a reasonable effort in field. However, the type of ecogeographical maps used as an input for the model must be chosen carefully according to the objective of the study. In our study the objective was to identify suitable areas in terms of geographic and geometric distribution. We therefore chose an unclassified Landsat scene and fitted this to the ENFA; the results from it seem to be satisfactory and robust (J. Sahlsten, unpubl. data). Furthermore, whether patches actually contain the species is not revealed by ENFA results, but the suitability of an area is probably valid (Hirzel et al. 2006). If the species is absent in a suitable patch, reasons other than those incorporated into the ENFA should be investigated.

The variation in patch size in our study is rather high, which may indicate a highly fragmented landscape. When a landscape is strongly heterogeneous, the matrix surrounding the patches should be considered in addition to the patches themselves in the analysis. The matrix has not been considered in this study because the matrix mostly affects dispersal. Since knowledge on cost of dispersal due to encounters with different landscape features was unknown, the matrix was considered homogenous. Furthermore, the study area is not isolated from its surroundings and it should be noted that patches at the periphery of the study area were probably affected by patches outside the study area. Thus, presence in peripheral patches is probably to some unknown extent explained by connection to patches not included in the model. This fact may lend a rescue effect to certain patches that could lead to a decreased extinction probability. However, due to the size of the study area, the distribution of patches and the ecology of hazel grouse, the impact of this unknown factor on the whole study area is small enough to be ignored.

Although, true absence data may be difficult to verify (MacKenzie et al. 2002), the high response frequency with the census method used in our study implies that there is no strong error or bias in presence/absence data. One could argue that the effort was different among patches due to differences in size. However, since the census was made by point-counting with a defined time period at each stop, the effort is the same per unit and different patch areas should not cause major bias to the result in either direction.

The α estimation may be problematic due to the difficulties in finding a value that is ecologically realistic. The mark-recapture methods used by Hanski (1994) can be used to estimate α, however, this method does not eliminate the problem of Euclidian distance (i.e. distance in a straight line). Therefore, friction maps or least cost distance would be a more accurate measure of distance if available (Ray et al. 2002). But knowledge of the cost to or mortality of hazel grouse when dispersing is not sufficient to create friction maps or to calculate cost distances. Hence, the Euclidian distance found by Montadert & Léonard (2006) stands as a good approximation of α for the model. It could be argued that the patches < 20 ha which were removed may function as stepping stones. However, our limited knowledge of dispersal behaviour did not allow us to incorporate stepping stones into the model. Hence, patches in our model were only included in terms of colonisation and extinction. It was assumed that the probability for colonisation was zero, and extinction probability was one when patch size was below a critical area coinciding with the minimum area for a home range, implying these small patches did not exist. Of course, it would have enhanced the model if small patches could have been included as facilitators of dispersal, and efforts to find good measures of dispersal distance should be a priority in further studies. In general, this seems to be lacking in many studies of species with large home ranges.

The area measurements used in our study are not claimed to be optimal, rather they indicate that the measure of area as an estimation of extinction probability should not be taken for granted. It should also be mentioned that the effect of the perimeter-area relationship probably has a threshold at the point where the area becomes large enough to sustain a population, regardles of the convolution of its edges. On the other hand, with such large patches in a system, the metapopulation approach is probably not appropriate and a source-sink or mainland-island system would be more suitable.