Active Management of Protected Areas Enhances Metapopulation Expansion Under Climate Change

There is a need to adapt biodiversity conservation to climate change, but few empirical studies are available to guide decision‐making. Existing networks of protected areas (PAs) have been preferentially colonized during species’ range expansions, but this could be due to their original habitat quality and/or to ongoing management activity. Here, we examine how PA status and active conservation management have influenced the range expansion of a butterfly Hesperia comma through fragmented landscapes. PAs under active conservation management were over three times more likely to be colonized than unprotected, unmanaged sites of the same basic vegetation type. Conservation action also increased the survival rate of existing populations inside and outside of PAs. We conclude that PAs facilitate range expansions by preventing habitat degradation and encouraging active conservation that improves habitat quality, and that conservation interventions on nondesignated sites also have a role to play in adapting conservation to climate change.


S1: Null model structure
We used the top-fitting models from a previous study (Lawson et al. 2012) as "null models" for the colonization and survival of Hesperia comma populations in this study. Terms indicating management and protection status were subsequently added to assess their effects over and above those explained by other variables. The structures of the null models are given below, and incorporate the following variables:  Patch area in hectares: the area of suitable habitat <10 cm in height was calculated for each patch by calculating the total area of each patch (using digitized patch polygons) and multiplying this figure by the proportion of the patch which contained turf less than 10 cm in height (estimated during field surveys).
 Host plant cover, representing the proportion of turf <10 cm in height which contained suitable Festuca ovina host plants.
 Bare ground cover, representing the proportion of turf <10 cm in height which contained bare ground. Bare ground is important for H. comma because it heats up more than longer vegetation in direct sunlight and provides warm microclimates for egg-laying (Davies et al. 2006). The survival model presented below includes a squared bare ground term to account for the fact that patches with very high bare quantities of bare ground are less suitable for H. comma populations (Lawson et al. 2012).
 Direct connectivity, representing the expected number of adult butterflies to immigrate into the focal patch from surrounding habitat patches.
 Indirect connectivity, representing the connectivity of the focal patch to empty but suitable habitats, which facilitate immigration events.
The two connectivity measures were calculated using the following formula: Where i is the focal patch and j all other patches, which have area A j and are separated from i by distance d ij . Here, A j is effective area <10 cm (ha) of patch j, and d ij edge-to-edge distances between patches i and j (km). α (a negative exponential dispersal kernel) and b (a scaling function for patch emigration) are estimated from a previous study (Wilson et al. 2010). For direct connectivity, p=1 for occupied patches and 0 for unoccupied patches. For indirect connectivity, p j is calculated is ∑ ≠ (i.e. the connectivity of each patch j weighted by its area A j ). The original analysis which demonstrates the importance of the variables in the above variables for the establishment and survival of H. comma populations can be found in Lawson et al. (2012). Colonization Where indicates the probability of colonizing patch , indicates the "direct connectivity" of patch to patches that were occupied in 2000, and indicates the "indirect connectivity" of patch I, reflecting the availability of suitable but uncolonised habitat surrounding the patch (see Lawson et al. 2012), and indicates the proportion of the patch that was covered by the host plant Festuca ovina. is a binary variable indicating whether patch was colonised between 2000 and 2009. Survival ( ) = + 1 + 2 + 3 + 4 + 5 + 6 2 ~( ) Where indicates the probability of survival in patch , indicates the "direct connectivity" of patch to patches that were occupied in 2000, indicates the areal extent of the patch, indicates the solar index of the patch (a combination of aspect and slope), indicates the "macroclimate" of the patch (mean daily August maximum temperature from 2000-2009), and indicates the proportion of the patch that was bare ground. is a binary variable indicating whether the population in patch remained present in 2009.

S2: Investigating the influence of spatial autocorrelation
Spatial autocorrelation can have important effects on the conclusions derived from models of species distribution dynamics (Dormann et al. 2007). We examined our data for spatial autocorrelation, and assessed the spatial scale at which spatial autocorrelation occurs, using a semivariance analysis (Meisel and Turner 1998). Because the analyses presented in this manuscript are based on null models (see section S1) which incorporate spatially-explicit biological processes (e.g. connectivity effects), we used a semivariance analysis to examine the spatial patterns in the residuals of these null models. The results of this analysis are presented in Figure S1.
In the null models for colonization ( Fig. S1a), the semivariance reaches a sill (leveling-off point; Meisel and Turner 1998) between 10 and 15 km, possibly extending up to 20km. This indicates that spatial autocorrelation is strongest within this range (most notably at distances within 10km; Fig  S1a). For the survival null models (Fig. S1b), a sill is reached within 10km, but there is an apparent spike in semivariance at just below 25km. Note, however, that our dataset contains few pairs of patches separated by this distance range (number of pairs = 91), so our estimate of the semivariance at this point is somewhat uncertain. Overall, we can be confident that spatial autocorrelation in the residuals of the null survival model occurs within 25km. For this reason, we choose to examine the effects of spatial autocorrelation on our conclusions regarding management and protection status at a distances ranging from 5 to 25 km (see below).
To assess the extent to which our results were robust to spatial autocorrelation effects, we repeated our analyses using generalized linear mixed models. Instead of a single intercept for all sites ( in the equations above), we fitted a random intercept which grouped patches within grid squares of sizes ranging from 5km to 25km, using the lme4 package in R (Bates et al. 2011; R Development Core Team 2011). To calculate patch groups, we constructed R code (R Development Core Team 2011) which superimposed a grid of each resolution (5km, 10km, 15km, 20km or 25km) across the British Isles, with the origin (southwest corner with spatial coordinates [0,0]) taken to be the origin of the UK Ordnance Survey National Grid. We then classified patches into grid square groups based on whether they were within the same grid square. As such, the fixed intercept was replaced with a random intercept: Where is the intercept for patches in grid cell group , is the mean probability of colonization or survival across all patches, and 2 represents the variance in the probability of colonization or survival among grid cells. Table S3 displays model selection tables for colonization models with and without spatial autocorrelation effects. The rank order of the best models and the direction of predictions remain consistent whether or not a random intercept is used, indicating that the findings of this analysis are robust to spatial autocorrelation. Table S4 displays model selection tables for colonization models with and without spatial autocorrelation effects. The finding that primary management improves the probability of survival remains consistent whether or not the effects of spatial autocorrelation are considered (Table S4). However, the evidence that populations were more likely to survive in protected patches weakens once spatial autocorrelation is accounted for, such that the effects of protection status on survival may have been exaggerated by spatial autocorrelation effects.

S3: Detail on method for predicting colonization and survival probabilities
Table 2 (main text) displays predicted colonization and survival probabilities for a patch in each of the management categories. Because our models also incorporated effects of other patch attributes (e.g. patch size and connectivity; Lawson et al. 2012), we needed to choose values for these variables to produce colonization and survival predictions for an "average" patch. We chose to use the mean values based on all patches used in each analysis (colonization and survival). The values of these variables are given in Tables S1 and S2. Note that only variables which entered models (see section S1) are given.

S4: Relationships between management categories and variables in null models
Some management categories were associated with improvements in colonization and survival probabilities (Fig. 1, main text), but had apparently little benefit once important environmental variables were controlled for (Table 1, main text; see section S1 of this supplementary material for environmental variables included in models). To further explore why this might be, we investigated whether management or protection designation was associated with the environmental variables included in our null models. In Figures S2 and S3 we plot the distributions of environmental variables for patches used in (a) the colonization analysis (Fig. S2) and (b) the survival analysis (Fig. S3). In the following paragraphs, we briefly discuss differences in environmental variables among protection and management categories.

Colonization data
Amongst sites that were unoccupied by H. comma in 2000, patches that were close to existing H. comma populations (i.e. patches that had higher direct connectivity) and were in more wellconnected networks of habitat (i.e. had higher indirect connectivity; see Lawson et al. 2012) tended to be under primary management by conservation bodies, rather than voluntarily managed under agri-environment schemes or unmanaged. Our models therefore suggest that primary management greatly improved colonization chances of patches over and above the benefits of their higher connectivity (Table 1a, main text).
Voluntarily managed sites tended to have higher indirect connectivity than unmanaged sites, which might explain why we found little evidence for positive effects of voluntary management on colonization once the effects of connectivity had been accounted for.
There was no overall tendency for protected sites to have higher connectivity (direct or indirect) or host plant cover than unprotected sites, supporting our conclusion that protected areas improved colonization independently of these variables (Table 1a, main text).

Survival data
Amongst habitat patches occupied by H. comma in 2000, voluntarily managed sites tended to be larger and in more well-connected habitat networks (higher indirect connectivity) than either managed or unmanaged sites. Thus, patch size and connectivity variables may have played a role in influencing land managers' decisions to "opt-in" to agri-environment schemes (AES), and could explain why we found no positive effect of AES once these variables had been accounted for (Table  1b, main text).
On average, sites with higher direct connectivity were more likely to be protected as Sites of Special Scientific Interest (SSSIs). This reflects the fact that the 2000 distribution of H. comma was concentrated around protected areas. Protected sites also tended to be found in warmer regions of Britain (i.e. with higher mean August maximum temperatures). Both of these variables may have exaggerated the impacts of protected areas on population survival between 2000 and 2009 (Fig. 1b, main text), explaining why we found only relatively weak evidence that protection enhanced population survival (Table 1b, main text) despite survival in protected areas being higher than in unprotected areas (Fig. 1b, main text).

Supplementary tables
Variable Value Direct connectivity 2.0 Indirect connectivity 4.9 Host plant cover (%) 16