Host plant density and patch isolation drive occupancy and abundance at a butterfly's northern range margin

Abstract Marginal populations are usually small, fragmented, and vulnerable to extinction, which makes them particularly interesting from a conservation point of view. They are also the starting point of range shifts that result from climate change, through a process involving colonization of newly suitable sites at the cool margin of species distributions. Hence, understanding the processes that drive demography and distribution at high‐latitude populations is essential to forecast the response of species to global changes. We investigated the relative importance of solar irradiance (as a proxy for microclimate), habitat quality, and connectivity on occupancy, abundance, and population stability at the northern range margin of the Oberthür's grizzled skipper butterfly Pyrgus armoricanus. For this purpose, butterfly abundance was surveyed in a habitat network consisting of 50 habitat patches over 12 years. We found that occupancy and abundance (average and variability) were mostly influenced by the density of host plants and the spatial isolation of patches, while solar irradiance and grazing frequency had only an effect on patch occupancy. Knowing that the distribution of host plants extends further north, we hypothesize that the actual variable limiting the northern distribution of P. armoricanus might be its dispersal capacity that prevents it from reaching more northern habitat patches. The persistence of this metapopulation in the face of global changes will thus be fundamentally linked to the maintenance of an efficient network of habitats.

high-latitude populations is essential to understand the factors that shape range limits and to forecast the response of species to climate change. Ultimately, this knowledge is crucial to inform the conservation of these vulnerable populations.
Beyond climate, several biotic and abiotic factors have the potential to drive population processes at latitudinal range margins and to determine the limits of species' ranges (Sexton, McIntyre, Angert, & Rice, 2009). For example, although there is ample evidence that many species are currently responding to climate change by shifting their distribution poleward (Chen, Hill, Ohlemuller, Roy, & Thomas, 2011;Gillings, Balmer, & Fuller, 2015;Hickling, Roy, Hill, Fox, & Thomas, 2006;Parmesan et al., 1999), the actual patterns of range shifts have been shown to result from a complex interaction between climate, biotic interactions (Van der Putten, Macel, & Visser, 2010), intrinsic species traits, and anthropogenic pressures (Jetz, Wilcove, & Dobson, 2007). In this regard, habitat fragmentation caused by human land use can be a key limiting factor. It may significantly impact range shift opportunities (Burrows et al., 2014) and can accelerate the extinction of isolated populations by disconnecting them from other suitable areas (Opdam & Wascher, 2004). Local habitat quality can also interact with climate variables in determining range boundaries and the response of populations to climate change (Kleijn et al., 2010;Nicolè, Dahlgren, Vivat, Till-Bottraud, & Ehrlén, 2011;Seabrook et al., 2014). Biotic interactors (prey, predators, or hosts) are important determinants of habitat quality and can strongly influence population performance and species distribution (Louthan, Doak, & Angert, 2015). They can be so important that a mismatched response to climate change can limit range shifts that would have otherwise occurred if species tracked solely their climate niche (Pelini et al., 2009;Schweiger, Settele, Kudrna, Klotz, & Kuhn, 2008).
One approach to determine the relative importance of these factors is to study their impact on the dynamics of high-latitude populations.
In this regard, population responses to different microclimates offer an indirect but useful assessment of the climate-related reaction norm of the species (Lawson, Bennie, Hodgson, Thomas, & Wilson, 2014;Thomas, 1993). If microclimate is the main driver of the abundance or distribution of these marginal populations, it would suggest that the species' range is likely limited by climatic factors. We can thus expect the species to react to climate change by shifting its range poleward (Bennie et al., 2013). If, instead, habitat fragmentation is already limiting current patterns, restoring landscape connectivity may be pivotal to the conservation of these populations. Otherwise, it is unlikely that the species will be able to cope with climate change and expand its range unless new habitats are created. Similarly, if population demography and distribution at leading-edge margins closely depend on such biotic factors, the response of the focal species to future changes will largely be driven by the response of its co-occurring species (Gilman, Urban, Tewksbury, Gilchrist, & Holt, 2010;Van der Putten et al., 2010).
Whether these populations can be the starting point of climate change tracking thus depends on their fine-scale drivers of distribution and dynamics and on the species intrinsic habitat requirements.
Here, we examined the effect of variation in solar irradiance (used as a proxy for potential microclimate), patch quality, or connectivity across a network of patches at a butterfly's northern range margin. We used a habitat network including a large part of the total population of Oberthür's grizzled skipper (Pyrgus armoricanus) ( Figure 1) in Sweden as a model to infer the importance of microclimate, patch quality, and connectivity on the regional distribution and local abundance at a species' northern range margin. The dynamics of P. armoricanus in southern Sweden has been proposed to be driven by metapopulation processes (Öckinger, 2006). Following metapopulation theory (Hanski, 1998), we expect the spatial configuration of patches-their area and their degree of isolation from the surrounding patches-to be an important driver of P. armoricanus abundance and habitat occupancy. This would highlight the importance of dispersal opportunities and thus the decisive impact of human land use on the current conservation status of the species and on its future response to climate change. Similarly, as vegetation structure is known to affect butterflies in general (Kruess & Tscharntke, 2002), and this species in particular (Eilers, Pettersson, & Öckinger, 2013), we also investigated the effect of variable grazing intensities on interpatches variation in occupancy and abundance. Moreover, we know from a previous study that microclimate influences the choice of oviposition sites in this species (Eilers et al., 2013). This factor is thus likely to affect the observed probability of habitat patches to be occupied and the population size they can sustain. Finally, owing to the fact that the presence of P. armoricanus in a grassland patch is generally closely linked to the availability of its host plant species (Öckinger, 2006), we also assessed the effect of the density of hosts. By ranking the importance of each of these factors in explaining variation in occupancy and abundance among habitat patches, we can gain insights into the population processes acting at high-latitude range margins and more specifically the potential response of this species to climate change.

| Study area and data collection
Our study species, Oberthür's grizzled skipper (P. armoricanus) (Figure 1), has a wide but fragmented distribution throughout North F I G U R E 1 Adult Oberthür's grizzled skipper (Pyrgus armoricanus). Photograph by Theresia Widhalm and Alexander Neubauer Africa and Europe. Its northernmost populations are located in southern Scandinavia (Sweden and Denmark), in a relative isolation from other populations in western and central Europe (Kudrna et al., 2011).
Its habitat consists of seminatural grasslands which host the specific plant species where it lays its eggs and on which its larvae feed.
Scandinavian populations are known to select primarily Filipendula vulgaris and Helianthemum nummularium (Christensen, 2000;Eilers et al., 2013), two species that are known to occur in a fragmented distribution up to 600 km further north (Hultén, 1971). The species has two generations per year: Spring generation adults fly from mid-or late May to mid-June, and the summer generation flies in August. In addition to Denmark, its Scandinavian distribution is restricted to a small area of ca. 30 × 20 km in southern Sweden where it occurs in a network of small and fragmented patches (Öckinger, 2007). There exist no records of historical occurrences further north (Eliasson, Ryrholm, Gärdenfors, Holmer, & Jilg, 2005;Nordström, Opheim, & Valle, 1955).
The habitat patches analyzed in this study (min area: 0.028 ha, max area: 14.71 ha), defined as patches of dry unfertilized grasslands with the presence of at least one of the host plants F. vulgaris and H. nummularium (Öckinger, 2006), were located in the core of this system, mainly around the town of Tomelilla (Figure 2). Adjacent habitat patches were defined as discrete if separated by at least 50 m of divergent vegetation (often arable fields or agriculturally improved grassland) as recommended by Ojanen, Nieminen, Meyke, Poyry, and Hanski (2013). Sites were monitored by slowly walking a transect (10 m width) covering the entire area of the patch. As transects were designed to allow the entire patch area to be surveyed, the transect length was proportional to patch area, but, the sampling effort per unit area was constant.
Counts should thus reflect butterfly abundance within patches and not only transect length. All observed adult P. armoricanus individuals were recorded. If necessary, species identity was confirmed by capturing individuals with a handheld net. As we observed that the abundance of the summer generation is on average almost three times higher than the spring generation (Table 1), we considered that processes acting in each of them might be different (see, for example, Roy & Thomas, 2003) and thus analyzed each generation separately. Therefore, we derived for each habitat patch and each generation: occupancy and the average and variability of abundance. Occupancy was expressed, separately for spring and summer generations, as the number of years a patch was occupied divided by the number of years it was surveyed in this generation. As imperfect detection can bias estimates of butterfly abundance and occupancy, it is sometimes advised to use multiple surveys per season to accurately estimate occupancy rates (MacKenzie F I G U R E 2 Location of the study area (shown as a cross in upper left inset) and spatial configuration of the 50 sites surveyed in this study et al., 2002). Instead, we chose here, due to time limitations, to maximize the number of patches visited each year rather than visiting each patch multiple times. However, we accounted for possible incomplete detection by considering a patch unoccupied only when no individuals were recorded during two consecutive surveys, whether these surveys occurred in the same year (spring and summer generations) or in different years (summer generation at year t and spring or summer generation-when only summer generations were surveyed-at year t + 1). This conservative estimate should ensure the robustness of our results even if some sites were mistakenly assumed to be unoccupied.
Abundance was measured only for patches that have been recorded as occupied in at least one survey, regardless of the generation. This means that, for example, if a patch was only occupied in one summer generation survey, it was still retained for the abundance analyses of the spring generations, where it thus had a mean abundance of zero.
Abundance variability was defined as the coefficient of variation of abundances across years.
We characterized each habitat patch by six variables that described habitat quality, solar irradiance, or spatial configuration of patches (Table 1). In each habitat patch, we estimated the density of and H. nummularium, although the former plant species is typically much more abundant (Eilers et al., 2013). Therefore, we pooled the cover of F. vulgaris and H. nummularium and used the averaged value over the 10 plots as a measure of host plant density per patch. All patches had a density of host plants between 0% and 15%, except one habitat patch that showed an exceptionally high host plant density of ca. 34%. We also categorized each patch by its grazing frequency according to three classes: sites that were never grazed during the whole period of survey (thereafter referred to as "never"), sites that were grazed every year ("always") and sites that were grazed or not depending on the year ("sometimes"). We used solar irradiance to characterize potential microclimate, a feature that is known to affect P. armoricanus oviposition site selection (Eilers et al., 2013). Solar irradiance depends on latitude, elevation, aspect, slope, and surrounding topography and reflects the level of energy that is received at a given point of Earth. We used a digital elevation model from the Swedish Lantmäteriet (2015), produced by laser scanning. The accuracy of the elevation model is 0.5 m, and the resolution of grid cells is 2 m.
Insolation was estimated as the total direct solar irradiance per 2 m grid cell per year (Wh/m 2 ), using the solar radiation function in the Spatial Analyst toolbox in ArcGIS 10.2 (ESRI Inc., Redlands, CA, USA), based on latitude, slope, aspect, and effects of shading from the elevations of surrounding cells. We calculated the mean and standard deviation of solar irradiance for each habitat patch to represent the average and variability in microclimate.
We also quantified for each patch the two main predictors of metapopulation dynamics: area and connectivity. The boundaries of habitat patches were defined in the field and digitized in a GIS. To account for the fact that transect lengths were dependent on patch area (see monitoring protocol above), we included area as a covariable in all analyses.
Patch areas were calculated using ArcGIS 10.2. As a measure of connectivity, we used the index S i developed by Hanski (1999) and calculated as S i = ∑ j≠i e −αd ij N j . S i estimates the connectivity at patch i, where d ij is the Euclidian distance between patches i and j (in meters) and N j the observed abundance at patch j (reflecting the number of potential emigrants from patch j), α being a constant describing how fast immigration probability from patch j decreases with increasing distance.
Here, α was set as 0.0034, corresponding to an average movement distance of 295 m, according to a previous analysis of mark-recapture data (E. Öckinger, unpublished data). When calculating connectivity, we took into account not only the 41 patches analyzed in this study, but all potential habitat patches within a 35 × 35 km square including the entire known distribution of P. armoricanus in Sweden (25 patches in addition to the 50 surveyed patches). As this index depends on abundance of all surrounding patches, we calculated S i for each year and generation, and used in further analyses its average value over all years for each generation separately. When abundance data were unavailable in a given year, we used the average abundance at this site across all years of survey. Similarly, for all sites that were not part of the annual monitoring, we used a rough estimation of abundance T A B L E 1 Average characteristics of the 50 sites surveyed. Response variables: mean (±SD) occupancy, average, and coefficient of variation (CV) of abundance per patch for each generation. Explanatory variables: mean (±SD) connectivity (depends on abundance of other patches, hence one value for each generation), area, host plant density, mean and standard deviation of solar irradiance per site, and number of sites in each category of grazing at these sites based on a systematic mapping of the species' entire distribution in Sweden in 2007 and 2010.

| Statistical analyses
We analyzed the effect of habitat quality, microclimate, area, and irradiance. In addition, we aimed to test whether the proximity of neighboring patches, and thus a higher immigration probability in a metapopulation model, could balance low habitat quality or unfavorable microclimatic conditions. Therefore, we also included in all models as explanatory variables the two-way interactions between connectivity and site quality variables (host plant density and grazing) or microclimate (mean and SD of solar irradiance), resulting in four additional predictors.
We adopted an information-theoretic approach (Burnham & Anderson, 2002) by computing models with all combinations of variables, and ranked them by their second-order Akaike information criterion (AIC c ). A multimodel inference was then performed by averaging all models whose cumulative Akaike weight was <0.95 (Burnham & Anderson, 2002). We extracted standardized averaged parameter estimates of all variables and interactions and estimated relative variable importance based on the sum of Akaike weights of all candidate models containing the variable. Multimodel inferences were run using the "MuMIn" package (Barton, 2013)

| RESULTS
Nine patches were occupied through all 12 years, nine patches were never occupied, and 32 patches were occupied during at least 1 year.
Patch occupancy thus ranged from 0 to 1 (mean = 0.51 ± 0.39 SD) when both generations where considered together, with roughly similar values in the spring (from 0 to 1, mean = 0.47 ± 0.39 SD) and the summer generations (from 0 to 1, mean = 0.52 ± 0.40 SD) (Appendix 1 and Table 1  For average abundance, results from model averaging showed highly similar responses for both generations (Table 2 and Figure 3b).
The most important variables explaining average abundance were patch area, connectivity, and host plant density (relative importance in all cases > 0.9). The interaction between connectivity and host plant density had also a high importance in the models (>0.8). As for occupancy, average abundance increased with patch area and host plant density or connectivity. Again, the interaction between the latter two revealed that highly connected patches were less sensitive to host den- T A B L E 2 Ninety-five percent model-averaged coefficients (±SE) and variable importance from linear models explaining for each generation (a) occupancy, (b) average, and (c) coefficient of variation of abundance per patch. For grazing regime, the "never" category is taken as reference. Variable importance > 0.5 is highlighted in bold font, and the three highest absolute coefficient value for each model is displayed in italic font For the variability in abundance among years, models led to some notably different responses between generations (Table 2). In the spring generation, the variability in abundance tended to be reduced in large well-connected patches with a high host plant density (relative importance always = 1) and regularly grazed (0.72), although again the effect of plant density was reduced in highly connected patches (relative importance for interaction term = 0.99, Figure 3c).
Model-averaged coefficients revealed that host plant density, alone (coefficient = −1.13) and in interaction with connectivity (coefficient = 1.27), had the highest effects on site variability during first generation. In contrast, the relative importance was more evenly

| DISCUSSION
Although limiting factors may change in space and time (Lawson, Bennie, Thomas, Hodgson, & Wilson, 2012), understanding the determinants of distribution and population dynamics at species range margins is fundamental for our ability to predict biodiversity responses to climate change. We demonstrated that the occupancy, abundance, and population variability at the northern range margin of the butterfly P. armoricanus are mainly driven by patch area, connectivity, and host plant density, while solar irradiance only had an impact on patch occupancy. Large and well-connected habitat patches tended to be more often occupied, to have on average larger and more stable populations, and to display a higher abundance. The strong and consistent effects of patch area and connectivity show that the spatial configuration of habitat is the most important factor for population persistence at the climatic range margin of this butterfly. This pattern is congruent with metapopulation theory which predicts that occupancy and turnover rate are driven by patch area and site isolation (Hanski, 1994(Hanski, , 1998.

F I G U R E 3 Partial relationships
showing model averaging predictions for (a) occupancy, (b) average, and (c) coefficient of variation of abundance across the range of host plant density observed in the patches. As models included the interaction between host plant density and connectivity, we plotted predictions for three levels of connectivity corresponding to 0.25, 0.5, and 0.75 quantiles (connectivity increases from lighter to darker colors). Plots are shown for models built with spring (red, left) and summer generation (blue, right) data, and observed values are displayed as triangles and circles, respectively Several studies have highlighted the importance of microclimatic conditions for invertebrate populations, especially for populations near the climatic limits of the species (e.g., Bennie et al., 2013;Turlure, Choutt, Baguette, & Van Dyck, 2009;Wilson, Davies, & Thomas, 2010). We found that patches that received a high solar irradiancepresumably reflecting a warmer microclimate-were more frequently occupied than less sunny habitat patches, confirming previous findings that microclimate is an important aspect of habitat quality for this species (Eilers et al., 2013) (Suggitt et al., 2011). A previous study found that the availability of host plants situated in a warm microclimate could predict local population sizes (Eilers et al., 2013). With warmer global temperatures, it is possible that this restriction to sites with a warm microclimate can be relaxed, and allow for colonization of a wider range of habitats at a regional scale, as has been observed for other butterflies at the climatic margins of their distributions (Pateman, Hill, Roy, Fox, & Thomas, 2012;Pateman, Thomas, Hayward, & Hill, 2016;Wilson et al., 2010). So far, however, this has not been observed for P. armoricanus in our study region.
The density of larval host plants, F. vulgaris and H. nummularium, appeared as a major driver of occupancy, abundance, and population variability. Larval host plants have long been known to play a vital role in the dynamics of butterfly populations (Hanski & Singer, 2001;Koh et al., 2004). As such, host plants' dynamics and distribution can strongly affect butterflies' response to climate change (Araújo & Luoto, 2007) and predicted range shifts by host plants could directly lead to changes or reductions in the distribution of herbivorous insects (Romo, Garcia-Barros, Marquez, Moreno, & Real, 2014;Romo, Silvestre, & Munguira, 2015). On the contrary, a mismatched response between a butterfly and its host may make it unable to track climate change (Pelini et al., 2009;Schweiger et al., 2008), even if host switching has also been documented (Pateman et al., 2012;Thomas et al., 2001 (Hultén, 1971). This leads us to believe that a potential expansion by the species to habitat further north is restricted by either microclimatic favorability or connectivity.
Connectivity had a particularly high importance in almost all models, especially in interaction with host plant density. Considering the effect of connectivity alone, our results reflect the fact that sites that are spatially isolated and surrounded by low-abundance sites have generally a lower and more variable abundance and are less often occupied. This suggests that the current distribution of P. armoricanus could be constrained by a too high isolation of habitat patches further north, hence preventing their colonization due to limited dispersal. The effect of the interaction between host plant density and connectivity is also interesting; it reveals that the positive effect of host plants density on abundance or occupancy-and to a lesser extent site stabilitytends to decrease as connectivity increases. It can likely be interpreted as a rescue effect (Brown & Kodric-Brown, 1977) which allows lowquality sites to be sustained by migration from surrounding patches and reduced population fluctuations in well-connected patches. This property of metapopulation dynamics favors the persistence of a local population by decreasing extinction probability or by supporting population size through regular immigration events, even if the local conditions are suboptimal (Gonzalez, Lawton, Gilbert, Blackburn, & Evans-Freke, 1998;Gotelli, 1991). Because patches are separated by unsuitable intensively managed fields, the isolation of habitats is in the present case strongly driven by agricultural practices which determine the connectivity of the habitat network and will in this respect be key to the long-term persistence of this P. armoricanus population.
Our results challenge the still-common view that latitudinal edges are limited purely by climatic factors (Pearson & Dawson, 2003;Woodward, 1990), but concurs with modern niche theory that assumes actual distributional limits to be formed by an interaction between abiotic factors (fundamental niche, mainly climate), biotic interactions (here larval host), and dispersal (Soberón, 2007;Soberon & Nakamura, 2009).
Better understanding of the actual position of this population relative to the species niche could be gained by investigating more closely its climatic tolerance, for example experimentally, or by comparing processes acting in various parts of its range (Lawson et al., 2012). In conclusion, it appears that the regional distribution and abundance of the northernmost population of P. armoricanus are mostly dependent on the availability of habitat patches with a high density of its larval host plants and on its capacity to disperse between such habitat patches. An effective conservation management strategy for this species should thus act both at the patch and landscape scales. First, habitat quality of already suitable patches must be maintained by ensuring the continuation of extensive grazing practices that provide an adequate vegetation structure and density of host plants. Second, the network of habitat patches must be kept dense enough to allow the long-term metapopulation viability.
More generally, the persistence of many species in the face of climate change will be fundamentally linked to the maintenance of an efficient network of habitats (Hodgson, Thomas, Wintle, & Moilanen, 2009). It can be achieved by preserving the existent connectivity between habitats, but also by creating or restoring habitats to facilitate range shifts across a fragmented landscape (Hodgson, Wallis, Krishna, Cornell, & Isaac, 2016;Hodgson et al., 2011). In addition, enhancing local habitat quality at climatic range margins has been shown to be an effective alternative strategy to facilitate species expansion under climate change, as it secures vulnerable marginal populations and increases the pool of potential migrants (Lawson et al., 2012). Conservation planning should also take into account the current and future properties of the landscape matrix to assist species in tracking their favorable

APPENDIX 1
Frequency histograms of occupancy, average abundance, and coefficient of variation of abundance, for the spring (red) and summer (blue) generations. Overlapping areas appear in purple.

APPENDIX 2
Models ranked by AIC c , explaining (a) occupancy, (b) average, and (c) coefficient of variation of abundance, for each generation. Only models used in model averaging are shown, that is, those whose cumulative weight was < 95%. Full models (including all variables and interactions) and null models (including only intercept) are also shown for comparison. Variables are abbreviated as follows: A: area; C: connectivity; HD: host plant density; G: grazing intensity; SRM: mean solar irradiance; SRSD: standard deviation of solar irradiance. (Continues)